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CHAPTER 1: INTRODUCTION 


The purpose of this thesis is to study the stability characteristics of a pair of coupled van 
der Pol equations. In order to fulfill this purpose, several different concepts must be 
discussed prior to actually looking at the behavior of the coupled equations. 


The first issue to be dealt with is the van der Pol equation itself; this topic will be covered 
in this chapter with details of the solution to the van der Pol equation via perturbation 
theory given in Chapter 2. 


Chapter 3 will explore the derivation, existence and computation of Padé approximants 
which will be used to model behaviors of the van der Pol equation and the coupled 
equations. If the reader is already comfortable with the characteristics of Padé 
approximants, Chapter 3 can be skipped. 


In Chapter 4, the model for the coupled van der Pol oscillators will be examined. The 
scope of the problem of characterizing the stability behaviors of the coupled equations will 
be examined. To that end, the Zero Mean Damping (ZMD) surface will be introduced and 
the variational equations to the coupled oscillators will be derived. Some implications of 
Floquet theory on the variational equations will also be explored. 


Chapter 5 will detail the solution process of the variational equation derived in Chapter 4. 
The stability transition curves along the ZMD surface will be introduced and the method 
for determining them will be presented. 

Chapter 6 will summarize the results of Chapter 5 and discuss some ongoing work and 


computational difficulties encountered in determining the stability transition curves. 
Further studies will also be discussed. 


THE VAN DER POL EQUATION 


The van der Pol equation: 


i-e(1—u’)u+u=0 (1-1) 





2 
was introduced by Bismarck van der Pol in a paper published in 1927 [29,30] to model the 
non-linear resistance of certain circuits including vacuum tubes. Since that time Equation 
(1-1) has become a popular mathematical model for limit cycle oscillators. 


A limit cycle oscillator, is an oscillator which exhibits self-sustained periodic behavior. It 
is non-conservative, due to a damping term, and non-linear [2,7]. In the case of the van 
der Pol equation, the non-linearity is also due to the damping term. Limit cycle behavior 1s 
present in many mechanical, electrical, and biological systems [3,4,5,6]. Stoker[21] and 
LaSalle[15] showed the existence of a unique limit cycle for Equation (1-1) for all e>0. 
The limit cycle in the phase plane {u, u’} is shown in Figure (1-1) for the case of s=0.9. 


The shape of the limit cycle for the van der Pol equation varies markedly with the value of 
the non-linear parameter, ¢. For e<<1 the limit cycle is well approximated by u=2 Cos(t), 
while for large & the limit cycle remains periodic but varies significantly from sinusoidal 
[10,28,29]. For large s, the van der Pol equation can be used to model relaxation 
oscillators such as a heart beating, a clock ticking or a fish swimming [6,11,22,28]. 


i) 


LC) 


Figure (1-1) Phase-plane for VDP 
Oscillator (¢=0.9). 


Applications of relaxation oscillators is the motivation for this thesis. If the behavior of a 
chain of linked van der Pol equations can be characterized, it may be possible to build an 
electro-mechanical system which uses coupled oscillators to produce control signals for 





3 
locomotion. As an example, consider a pair van der Pol oscillators in the relaxation limit. 
Since u(t) is characterized by a rapid transition between two relatively flat states, using the 
signal to drive a pair of servomechanisms may be a very inexpensive way to drive a bipedal 
robot. By coupling the oscillators, it may be possible to produce a smooth, responsive 
gait with little or no computer processing. 


Another possible application for coupled van der Pol oscillators could be modeling the 
synapses in a fish’s spine while swimming [6]. Here the chain of oscillators may be 30 or 
even 50 units long and modeling may be extremely complicated. 


In the first step to understanding the behavior of such models, it is desirable to find the 
stability characteristics of the coupled equations. In this thesis, the scope will be limited to 
only two oscillators which could be modeled as: 


x —e(1— x? )k +x =eA(y—x)+eB(y— x) 


1-2 

y—e(1—y’ y+ y =eA(x-y)+eB(x-y) a 

where €A represents the displacement coupling and €B represents the velocity coupling 

between the two oscillators. This model was proposed by Rand and Holmes [18] and 

further studied by Storti and Reinhall[19,20,22,27]. In future, examining the stability of a 

chain of more than two oscillators could be done by extending the techniques of Storti and 
Reinhall and this thesis. 


CHARACTERIZING THE VAN DER POL EQUATION 


The first step in understanding the stability of Equations (1-2) 1s to develop the solution to 
the van der Pol equation. Equation (1-1) has been studied extensively [1, 2, 7, 8] and a 
simple approach to finding solutions using Lindstedt’s method [14] is presented in Chapter 
Zz 
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CHAPTER 2: SOLUTION TO THE VAN DER POL EQUATION 


Perturbation theory is one of the most powerful techniques used in solving non-linear 
differential equations in vibration theory [14,17]. In concept, perturbation theory is very 
simple and can be explored by examining the van der Pol equation in time Tt: 


ii—e(1—u’ )a+u=0, u(0) = 0 (2-1). 


Using Lindstedt’s method, rescale time and let t= @ t. Changing variables, Equation (2- 
1) becomes: 


u”—e(l—u’)ou’+@’u=0, ~—u’(0)=0 (2-2) 


where the primes represent differentiation with respect to the rescaled time t. Now apply 
perturbation theory, by considering Equation (2-2) as a small perturbation of the linear 
system: 


u”’+u=0, u'(0)=0 (2-3) 


whose behavior is well understood. It 1s reasonable to assume that u(t) and @ may be 
represented by power series in €: 


u(t) = u,(t)+eu,(t)+e-u,(t)+- 
@=1+€@,+8'@, +: (2-4) 
Substituting Equations (2-4) into Equation (2-2) yields: 
uvt+eult+euyt--+(1+eo0,+e8'o, +:)’(u,+eu, teu, +::)= 
e{1—(u, +eu,+eu, +) dteo, +£°@, ++:)(u, te ul teu, ++) 


(2-5) 


Collecting like powers of ¢, the first few equations become: 





6) ica ae 


l, " ne 2 
E> -U, pu =u, —u, j= Zen 


ae Z 
ce. uv tu, =@,us(1—u,’) +us(l-u,’) — 2u,u,u, —2@,u, —20,u, -@,"U, 


(2-6). 


These equations are then solved, in order, where the right hand side is completely fixed by 
the solutions to the previous equations and determining the @;’s and any undetermined 
coefficients in the earlier u ,(t)’s to suppress secular terms. 


Secular terms are terms that if not canceled in some manner, would give solutions to the 
differential equation that grow with time. As an example, consider the linear Equation (2- 
4) with an added driving force: 


u” +u = Cos(t) (2-7). 


Since the driving force matches the natural frequency, u(t) will grow over time and 
become unbounded; the solution to Equation (2-7) is u(t)=1/2 t sin(t), which clearly grows 
over time. In applying perturbation theory, the goal is to select the undetermined 
coefficients in Equation (2-4) to eliminate terms which match the frequency of Equation 
(2-3). In particular, for the van der Pol equation, the coefficients of any Sin(t) or Cos(t) 
term must be suppressed when they appear on the right hand sides of Equations (2-6). 


As a demonstration, consider the first few orders of & in Equations (2-6). Starting with 
zeroeth order, the solution must be of the form: 


u,(t) =A, Cos(t) +B, Sin(t) (2-8). 


Now applying the initial condition in Equation (7), Bo=0 and Apo will be determined while 
solving the first order equation. Substituting x(t) into the first order equation and 
employing trigonometric identities: 


w’+u, = (A,Cos(t))'(1- (A,Cos(t))”) — 2 ,(A,Cos(t)) 
(2-9). 


| l 
= (- A, + 1 A3) Sin(t) 2 qh Sin(3t) ag 2 @ ,A,Cos(t) 


In order to suppress the secular terms, the coefficients of Sin(t) and Cos(t) must be set 
equal to zero; there are three possible values for Ay to do this, Ag=-2, 0 or 2 and @; must 
be zero. Selecting Ap=0 will result in the trivial solution for the zeroeth order equation, 





6 
but either 2 or -2 would be valid and give solutions to Equation (2-1) with opposite signs. 
For convenience, take the solution Ap=2. The right hand side of Equation (2-9)is now 2 
Sin(3t) and the solution is: 


] 
u,(t) = B,Sin(t) — q sin) (2-10) 
and imposing the initial condition of Equation (2-1) again: 
a, oe 
u, (t) = 7 sintt) ~ qoinGo . (22m) 


It serves little point to continue this procedure here, but it is worth noting that both the 
first order term of the solution, u(t), and the first order term in the frequency expansion, 
@i, were determined from the first order equation of Equations (2-6). This pattern will 
hold while solving Equations (2-6) at each order of &. Following this procedure, it is 
possible to find solutions for all u,(t)and ; up to a high value of 1; such a solution 1s 
considered an i” order solution to the differential equation [14]. 


This can be a tedious solution procedure, but the main advantage is that any degree of 
accuracy can be achieved. In particular, if the solution is needed to some high order n, 
this method will allow determining the solution of Equation (2-1) through order n with 
sufficient time and effort. Fortunately, this method (Lindstedt’s method) can be realized 
using a programmable, symbolic algebra package. The results outlined here were 
determined using a program written in Mathematica™ [9, 32]; the code to generate these 
results is found in Appendix A. This code takes advantage of some additional features of 
the van der Pol equation that will be discussed in the next section. 


SYMMETRY IN THE VAN DER POL EQUATION 


Examining Equation (2-6) with symmetry of the solution in mind, it can be seen that the 
right hand sides of the equations alternate in symmetry. In particular, the mght hand side 
of the zeroeth order equation is 0 which is even, therefore uo(t) must be even, which it 1s. 
The right hand side of the first order equation, considering the fact that u(t) 1s even, must 
be odd resulting in an odd solution for u,(t). This pattern continues and allows the rapid 
solving of each order of the power series solution. 


One additional property of the equation may be used to accelerate the solution of 
Equations (2-6). Since the left hand side of every order in € is the same, namely: 


LHS = u!’+u, (2-12) 
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then the homogeneous solution is always of the form: 
u,,(t) =A; Cos(t) + B, Sin(t) (2213) 
and due to the argument of the above paragraph, A;=0 if1 is odd or B;=0 if 1 is even. 
Even more fortunately, the particular solution for each order in ¢ is also predictable. The 


right hand side of Equation(2-6), after suppressing secularity, will always take the form 
of: 





(2i+1),A=2 
RHS = a, + an Cos(k t) V1 eEven 
k=3 
(21+1),A=2 ary. 
RHS= > b, Sin(k t) Vi eOdd 
k=3 
This results in particular solutions of the form: 
(21+1),A=2 a 
eA) = 3 2 ; + Cos(k t) V1 € Even 
~ (2-15). 


(2i+)A=2 pH 
x(0= 2) cytsSin(kt) Vi «Odd 


k=3 





The overall advantage of being able to write down these solutions, is that Lindstedt’s 
method may now be used without actually solving any differential equations. This 
significantly speeds up the solution generation for the van der Pol equation. 


After removing the need for using a differential equation solver, the next most time- 
intensive process left in solving each order of the van der Pol equation is expanding the 
right hand side of each equation and using the trigonometric identities to reduce them to 
the form of Equation (2-14). One very simple way to reduce the computational effort is to 
use Euler’s identity: 


e)' = Cos(t) + j Sin(t), or 


Bierce ae ee (2-16) 
Cost) =<—*— and Sin(t) = —<=— 
J 





8 
where j=vV-1. Further, in using Mathematica™, the expression Exp(j t) is very 


complicated, so it is advantageous to define another symbol, say Z = Exp(j t), and then 
replace each sine and cosine term on the right hand side of Equations (2-6) with: 


k -k k _o7-k 


ZL +z Z 
Cos(k t) a and Sin(k t) = oe (2-17) 
Using these expressions, the right hand side equations can be expanded much more rapidly 
and then regrouped back into terms of sines and cosines to implement Equation (2-15) and 
suppress secularity. 


The first few terms of the solution of Equation (2-1) are presented here: 


u (t) = 2 Cos (t) + 

3Sin(t) 1. 2, Cos (t) 3 5 

( "4 - 7 Sin (3b) +¢ (- 3 + <= Cos (3 t) - <= Cos (5 t)) + 
af Se 2 te a Sauer (a ‘ 
(- + 5 Sin (3 t)- = Sin (St) + —— Sin (7) +0.) 

















256 256 (2-18) 
with the frequency expansion given by: 
e 17e4 356 6788998 28160413 «1° 6 
w=t1-—t + t+ Ht ) 
16 ©3072» 884736 §=—- 5096079360 ~=——- 2293235712000 (2219) 


Using the short-cuts listed above, the code found in Appendix A was developed. The first 
ten terms in the solution of the van der Pol equation are found in Appendix C along with 
numerical solutions up through order 25. The solution to the van der Pol equation was 
determined analytically (no numerical rounding) up through order 75 using this code. 





CHAPTER 3: BRIEF DESCRIPTION OF PADE APPROXIMANTS 


WHAT IS A PADE APPROXIMANT? 


A Padé approximant is a rational function which approximates a power series [12,13]. It 
is represented by [L/P] where the numerator is of degree L and the denominator 1s of 
degree P. The Padé approximant [L/P] can be constructed from any power series of 
degree N: 


0 


f(z) = >.¢,z' =c, +¢,2+0,27 +--+, 2% + O(28"") (3-1) 
i=0 


where L + P <N, and takes the form: 


Apa eae 7: tae 
L/ P| =——2_>— + 0(2" 32 
Lie b, +b,z+b,z +--+b,z (2h) G2) 


and the Padé approximant’s Maclaurin series expansion must match Equation (3-1) 
exactly up through order L + P. In order to make [L/P] unique (since multiplying the 
numerator and denominator by any value would create another approximant with the same 
Maclaurin series), it is common to arbitrarily set bo=1 (or divide both the numerator and 
denominator by bo), thus making Equation (3-2): 


2 L 
dy tly Zed 2, heed 


ee 
tee 1+b,z+b,z7+--+b,z’ 


+O(z*?*") (3-3) 


Note: The power series of Equation (3-1) is assumed to be complex- 
valued and could represent a Maclaurin series expansion of some function 
f(z); however, the power series need not be complex-valued in order to 
work with Padé approximants and Equation (3-1) could just as easily 
represent the Taylor’s series expansion of a real valued function, or an 
experimentally obtained polynomial curve fit, or a power series solution 
obtained with perturbation theory. 





WHAT ARE PADE APPROXIMANTS FOR? 
There are three main reasons for constructing Padé approximants. 


The first reason to construct [L/P] is to produce a function which has a larger radius of 
convergence than the original power series. In some cases [L/P] may have a significantly 
larger radius of convergence than the original power series, or possibly even converge for 
all values of z. Even if the approximant does not have a larger radius of convergence, the 
approximant may yield information about how the function behaves for large values of z. 


This will be more readily seen in an example. 
EXAMPLE 3-1 
Consider a power series: 
f(x)=1-x+x?—-x?4-- 


which is the Taylor’s series for 1/(1 + x) has a radius of convergence R<1. 
But the an approximant for f(x) is given by 


l 
(es 


which is the original function and is finite for all x except x=-land 
converges for R<oo. One can see the behavior of f(x) and [1/1] easily in a 
graph (f(x) is the dotted line, [1/1] is solid) where f(x) includes the first 8 
terms in the series: 


E(x) & [iy 
/ 






| 

1\_ Taylor’s Series 
| 
| 


[1/1] Approximant 






0 Xx 
i 2 3 


Figure (3-1) Taylor’s series and Padé approximant for 1/(1+x). 





1] 
Here, the approximant obviously does a better job of estimating the 
function for x>1, since the approximate is exactly the same as the original 
function. 


The second reason for using Padé approximants is to speed up the evaluation of a given 
power series when working inside the radius of convergence. For the preceding example, 
it is less computationally expensive to find the value of 1+ x and its reciprocal than it is to 
evaluate the Taylor’s series. In a more general sense, comparing the evaluation of two 
polynomials of degree N/2 and taking their ratio is a less costly evaluation than calculating 
the value of an N degree polynomial, due to the expense of calculating higher powers in 
the N degree polynomial. 


The third reason for employing Padé approximants, which is closely related to the first, is 
to more accurately model the asymptotic behavior of a function. One of the worst 
features about a power series solution is that the highest order term will always dominate 
when x is large, so the series behaves as x‘. This may be unfortunate, especially if there is 
some indication of how the actual solution should behave. For example, if the function 
being modeled is known to decay as 1/x or is known to approach a constant value, then a 
Padé approximant should be more successful in predicting the asymptotic behavior by 
careful selection of L and P when constructing the [L/P] approximant. In the case of a 
function which should decay as 1/x, selecting P=L+1 will probably result in an 
approximant that also behaves as 1/x in the limit as x becomes large since the numerator 
will behave as x and the denominator as x"! In the case where the function should 
approach a constant, picking L=P will most likely produce an approximant with the 
desired behavior. (Example 3 below will explore some of this more clearly.) 


CALCULATING PADE APPROXIMANTS 


PROCEDURE 


Calculating a Padé approximant is a straight-forward process [13]. Starting by setting 
Equation (3-1) equal to Equation (3-3): 


2 L 
andwanzeiaa 2 4°°°ta., Z 
ee OO 2 eZ 


Next, multiply both sides by the denominator of the right side: 


(c, +C,Z2+C, Z +40, 2)\(1+b,z+b,z?+--+b,z") = P ; 
2 L 
a), +a,Z2t4,Z +°°°t+a,Z 





Multiplying Equation (3-5) out and collecting like powers of z, yields: 


Crp) Cy-po2 7" Cr Ds Cray 
Cio Cpe ee ee ae (3-6) 
Cr Cyepay 7°" Cap b, Crip 


where c; = 0 if }<0 for consistency, and 
Ay = Co 
ace ce 
ayo Cu aabRe. Bom 


fin A) 


a,=c, + DbjC,, 


i=1 


The linear system in Equation (3-6) may be solved for the b;’s by any suitable method, and 
then the b;’s used in Equations (3-7) to find the a;’s. Following this procedure, any order 
Padé approximant can be constructed for L + P<N. 


Equations (3-6) and (3-7) can be easily coded into any programming language or written 
as functions in many mathematics packages. 


EXAMPLES 


EXAMPLE 3-2 


Find [4/4] for the Taylor’s series approximation to Sin(x). Starting with 
the Taylor’s series: 


Rei) Xe : 
f(x)=x-F +a, + Olx ) 
and [4/4] is 
a, ta,xta,x’> +a,x +a,x" 
ee 


i+bx+b,x’? +0, x° + bee 
] 22 3 4 





Using Equation (3-6): 


l 0 -1/3! 0 fb, 1/5! ] [b,] [11/5880 
0 -1/3! 0 1/5! |b, 0 b, 0 
1/3! 0 1/5! 0 |b, ]~ J-1/7!}7]}b, |>} 3749 
OWS) oe iy 7 0 b, 0 
a, =0 
a,=1+0=1 


a, =0+0+0=0 
a, = —-1/64+0+3/494+0=-31/294 
a,=0+0+07070—0 


so the final solution is: 


ee 
me 
[4/4]=— a + O(x?). 
ae: 4 
1+ 749% + 5gg0* 


Note: Padé approximants do not perform much better than regular power 
series in predicting periodic functions; this example was included for 
numerical interest as will be discussed later. 


EXAMPLE 3-3 


Another example of finding a Padé approximant demonstrates how 
asymptotic behavior for Padé approximants may be significantly better than 
a normal power series. Consider the function: 





l 
l Ls 
f(x) = === with Padé approximant [3 / a nn 
1+x’ [22 x? a 
8 


Selecting a [3/4] approximant increases the likelihood of correctly 
modeling the 1/x behavior of f(x) for large x. Attempting to calculate the 
[3/4] approximant yields a [2/4] approximant which shows that not all 
degree approximants are necessarily unique. Below is a graph of f(x) (solid 





line), [3/4] (dotted line) and the eighth order Taylor’s series for f(x) (dash- 
dot). In this case, the [3/4] approximant does a significantly better job of 
modeling the original equation and tends to zero (as 1/x’) for large x. 


Z 


[3/4] Approximant 


| 
\ 
a: \— Taylor’s Series 


Figure (3-2) Taylor’s series and a [3/4] approximant to f(x) = 





STABILITY OF PADE APPROXIMANTS 


CONDITION NUMBER CONSIDERATIONS 


Following the procedure set out here for finding a Padé approximant includes solving the 
linear system represented by Equation (3-7) and Equations (3-7). While Equations (3-7) 
present no difficulties computationally, the solution to Equation (3-6) may be subject to 
significant loss of numerical accuracy. This is of particular concern when combined with 
the fact that the accuracy of Padé approximants is highly sensitive to the differences 
between coefficients in the approximant. (The need for very high numerical accuracy is 
discussed in detail in Reference [13], section 2.1.) 


To illustrate how significant the numerical results may be affected when solving Equation 
(3-6), one can examine the condition number (or a reasonable approximation for the 
condition number) of the matrix C: 


7 ie - 
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cL == ~*~ =? 
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Cy pe hn Cy 
E e : ne 
L-P+2 L—P+3 L+1 
C= |e ae (3-8). 
Cy Cr Pi ee Cisp-1 


To gain some measure of the magnitude of the condition number of this matrix, consider 
the calculations carried out in Example 3-2. Below is a table which shows how the 
condition number of C varies with the degree of the Padé approximant used to estimate 
the function f(x)=Sin(x) and the number of lost digits in accuracy is assumed to be 
Logjo(condition number): 


L=P Cond# Lost Digits 
4 4.72 10° 4 
5 1.66 10° 5 
6 AQT A: 8 
qi 3.40 10° 10 
8 1.74 10” 12 
9 2.46 10" 14 


10 2.08 10!” 17 


Table (3-1) Condition number and loss of 
accuracy for Example (3-2). 


The importance of finding an accurate method to solve Equation (3-6) 1s apparent with 
respect to these results, and it would be wise to check the condition number of C when 
employing Padé approximants. 


There is one obvious way to avoid the numerical problems discussed here. By employing 
one of the powerful math packages available which can do exact arithmetic such as 
Mathematica™ or Maple™, Equation (3-6) can be solved with no loss of accuracy when 
the coefficients c; are known as precise fractions. While it may seem optimistic to want to 
know the c;’s exactly, in many applications where Padé approximants are used, the power 
series being replaced is generated exactly from a technique such as perturbation theory. In 
particular, in non-linear vibration analysis, the solution to a differential equation may be a 
power series in the non-linear parameter, €, known to 50 or even 100 terms which are 
calculated with fractions and known precisely. This will be discussed more fully when 
examining applications of Padé approximants to vibration analysis. 
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What is a defect? 


In simple terms, a defect is a failure of a Padé approximant to accurately approximate a 
power Series in the neighborhood of a root of the denominator of the approximant. Or in 
mathematical terms, a defect occurs around x=xp when 1+b,xp+b,xpt+--+bpx5, =0. 
Borrowing on the language of complex analysis, the defect occurs where [L/P] has a pole 
at Xp OF Zp. 


How to handle a defect 


The presence of the pole in [L/P] may or may not be a result of numerical inaccuracy in 
calculating the approximant, and on a first look it seems reasonable that there would be 
cases where we would expect the poles to be present. The most obvious reason for the 
[L/P] approximant to have a pole would be if the original function, f(z), also has a pole. 
In this case, the behavior of the approximant will depend on the nature of the pole in f(z); 
if the pole is an essential pole (there is no finite n such that (z-zp)” f(z) is analytic at zp) 
then the approximant will always have a defect at zp, regardless of the degree of [L/P]. 
However, if the pole in f(z) is pole of finite order, it may be possible to remove the defect 
in the approximant by constructing [L/P] with L and P large enough. There is no 
guarantee that such an approximant is possible to construct, but taking higher degree 
approximants frequently removes defects when applied to perturbation theory. 


A defect may also occur in an approximant when f(z) has no poles. In this case, the defect 
could occur due to computational errors which can be significant as discussed in the last 
section. To handle defects of this nature, there are two basic approaches. One approach 
is to use a more accurate method to compute the solution of Equation (3-6). This may 
entail finding the coefficients in the original power series to greater accuracy, using a 
machine capable of higher accuracy or resorting to rational approximants of the 
coefficients of the original power series and then solving Equation (3-6) exactly. The 
second method would be to try calculating the approximant to /ower degree in the hopes 
that the C matrix condition number will be enough smaller to prevent the resulting 
approximant from having a defect. 


Another cause for a defect to occur is trying to approximate a complicated function by too 
low an order Padé approximant. If this is the case, merely finding a higher order 
approximant will remove the defect with little effort. Coupling this comment with the last 
point in the previous paragraph, it seems that trying various degree approximants (higher 
and lower) until finding one that is defect-free is a reasonable approach; however, 
checking the condition number of C is an intelligent starting point to determine if the 
difficulty is computational. 


More information on making the calculation of Padé approximants numerically stable can 
be found in Reference [13]. 
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APPLYING PADE APPROXIMANTS TO PERTURBATION RESULTS 


Having briefly explored one type of perturbation theory application, there are some 
obvious opportunities to apply Padé approximants to good ends. 


One of the most direct applications for Padé approximants in vibrations is to predict the 
frequency response of a system in terms of the non-linearity parameter . There are many 
different methods for estimating the natural frequency of slightly non-linear systems (¢ << 
1), but many are inaccurate if the system is tuned such that ¢ is slightly larger. Given that, 
it is reasonable to calculate a Padé approximant for the frequency of a non-linear system 
obtained with a perturbation technique in an effort to accurately predict the response 
frequency for larger s. As a starting point, it is assured that a straight power series for 
will fail to accurately predict the frequency an oscillatory solution as € increases except in 
very special circumstances. In the general case the frequency will not grow as s" which is 
the asymptotic behavior of the power series for large ¢. This can be seen by examining the 
power series expansion for the frequency of the van der Pol equation. 


EXAMPLE 3-4 


Returning to the van der Pol oscillator, Equation (1-1), and using Lindstedt’s method, the 
frequency can be found to high order in s. (This result was derived in Chapter 2.) The 
first few terms are: 


] 17 
@ =1-—>8* +> 8 + 


35. «678899 
16° | 3072~ 


884736. 5096079360 ~ 





°+O(e'’). 


Taking the [4/4] approximant: 


492389 , 433361 , 


+e, + E 
2698560 34541568 
661049 , 1285087 


+. —_—__—_—_—. er 
2698560" | 57569280° 


[4/4] = 


4 


Graphing both the power series for @ (dashed) and [4/4] (solid): 
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Figure (3-3) Power series and Padé approximant to 
the frequency of a van der Pol limit cycle. 


From experimental data, the van der Pol oscillator is in fact a softening system, meaning 
the frequency of the response decreases with increased non-linear effects. Also, just from 
examining the graph, the behavior of the [4/4] approximant seems to continue the pattern 
of the power series expansion, before ¢ became large enough to make @ start to grow 
quickly. Both point to the fact that the Padé approximant does a good job of extending 
the approximation for frequency response as € increases. 


LARGE & CONCERNS 


Another area where Padé approximants are very useful in vibration theory is in application 
to relaxation oscillators. Relaxation oscillators have numerous applications in biological 
systems including how the heart beats, how fish swim (modeling the synapses in the fish’s 
central nervous system), and have possible application as control circuits for robotics. The 
relaxation oscillator is modeled by a non-linear differential equation where the non- 
linearity is allowed to become large compared to the linear components in the system, 
meaning in the neighborhood of 10 or 100 times larger. Consider the relaxation limit for 
the van der Pol oscillator, Equation (1-1), which behaves as a relaxation for large &. As 
seen in Example 3-4, the Padé approximant for the frequency more reasonably predicts 
behavior in the relaxation limit. 


APPLICATION TO COUPLED VAN DER POL OSCILLATORS 


The main reason for presenting this information on Padé approximants is in applying the 
theory to the stability transition curves that will be derived in Chapter 5. For the 
remainder of this thesis, some basic familiarity with Padé approximants will be assumed as 
they will be used to estimate the frequency of the limit cycle of Equation (1-1), the Zero 
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Mean Damping Surface described in Chapter 4, and the Stability Transition Curves found 
in Chapter 5. 





CHAPTER 4: COUPLED OSCILLATOR ANALYSIS 


PROPOSED MODEL FOR COUPLED OSCILLATORS 


As discussed in the introduction, the pair of van der Pol oscillators (x and y) will be 
modeled by two van der Pol equations coupled with a linear spring and a linear damper. 
The coupling will be “weak” in the respect that the spring and damper will both be order 
€; in particular the stiffness constant will be ¢A and the damping constant will be €B where 
€ is the same non-linear parameter as in the original oscillator. The coupled equations take 
the following form [18]: 


¥ —e(1-x?)k+x =eA(y—x)+eBly—x) 


(4-1) 
¥—e(1-y’ yt y =eA(x-y)+eBlx-y) 


where dots represent derivatives with respect to time T. 


The object is to determine the existence and stability of solutions to Equation (4-1). In 
particular, it is desired to find what values for eA and ¢B result in orbitally stable 
solutions. 


In this context, orbitally stable is most easily described in the phase plane. Consider the 
phase plane of the original van der Pol equation plotted in Figure (1-1); the limit cycle 
presented is orbitally stable. In mathematical parlance, for any small perturbation of the 
system away from the described limit cycle, there exists a finite distance 5 such that the 
maximum distance between any point on the disturbed path and the original limit cycle 1s 
less than or equal to 6 [14]. Understanding the concept of orbital stability is essential to 
comprehending the behavior of coupled oscillators. 


By inspection, an exact in-phase motion of the coupled oscillators 1s possible with 
x(t)=y(t). In that case, the coupling terms identically vanish and x(t) and y(t) will have 
to satisfy the original van der Pol equation as found in Chapter 1 with x(t) = y(t) = u(t). 
In this situation, the behavior of each of the coupled oscillators 1s exactly the same as the 
original solution, namely each oscillator would exhibit an orbitally stable limit cycle in its 
own two dimensional phase plane, as seen in Figure (1-1). 
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In the same vein, it is possible to take advantage of the known properties of u(t) by 
considering x(t) and y(t) as functions which vary from u(t) by small amounts. Treating x 
and y as variations away from the solution of the van der Pol equation clarifies the 
meaning of “stability” as applied to the coupled equations of Equation (4-1). While 
studying the variations away from u(t), determining the stability of x and y is equivalent to 
studying the origin in the variational space. This will be addressed in more detail after 
deriving the variational equations. 


DERIVATION OF THE VARIATIONAL EQUATIONS 


Treating x(t) and y(t) as small variations from u(t), the stability of the solutions of 
Equation (4-1) can be determined by characterizing the behavior of the variations. 
Consider both functions as small variations from the limit cycle u(t) [28]: 


x= => x=€+u 
Ve => v= ew 


$= 
ate (4-2) 


where € 1s the variation of x(t) from u(t) and 7 is the variation of y(t) from u(t). 
Replacing x and y in Equation (4-1) and rearranging terms yields: 


é — e(1-&? - 2ue—u2)E + (1+ 2eun)é + {i e(1-u?)u+ ul = 
eA(n-&) +eB(n-8) 

-e(1 =n? - 2un—u’)y + (1+ 2eud)y + {ii - e(1-u?)a +u} = 
eA(é — n) os sB(é _ i) 


(4-3) 


where the terms in braces are the original van der Pol equation and are identically equal to 
zero, resulting in: 
é ~ e(1 —€* —-2ué - u’)é +(14+ 2euuJé = sA(n ~ 4) + eB — é) 


4-4). 
ij e(1— n? — 2un-u?) + (1 +2eua)n = eA(E - n) +eB(E - A) Pa 


In the interest of studying only small variations from the limit cycle, non-linear terms in € 
and 1 are neglected: 


Oe 
ea = © a= -——= ‘ 
= eS 

ee => ap * 
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E — e(1 — u’ Je +(142eun)é = eA(n ~ é) + sB(i - é) 


. . e). 

n- (1 ~ u’)q +(1+2euu)n = eA(é - n) +eBé = i) 
In order to simplify the coupling terms appearing on the right hand sides of Equation (4- 
5), two new equations can be formed by adding and subtracting the two Equations (4-5) 
and letting: 


Wee ae 
oa (4-6). 
v= c= 1 
So Equations (4-5) become: 
Ww —e(1-u?)w+(1+2eud)w = 0 
(4-7). 


¥—e(1-2B-u2)¥ + (1+ 2euu + 2eA)v = 0 


By inspecting Equations (4-7) it can be seen that both equations are exact differentials 
which can be integrated with respect to time to yield: 


w—e(l-w)w+we=k, 


(4-8). 
v —e(1-2B-u?)v +(1+2eA)v=k, 


While studying stability, the solutions for v(t) and w(t) can be recast to further simplify 
Equation (4-6). The particular solutions to Equations (4-8) are constants, in particular: 


ax, (4-9). 
ve 1426A 





These solutions do not vary with time and can be removed by a simple coordinate 
transform in the variational space: w> w-w, and v-v-v,. Such a transformation will not 
alter the stability of the solutions, since stability is controlled by the exponential behavior 
of the functions. Then the varational equations become: 


w—e(1-u’)w+w=0 


(4-10). 
v —e(1-2B-u’)v+(1+2eA)v = 0 
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In order to study the stability of the coupled van der Pol oscillators, both of Equations (4- 
10) must be explored and characterized. In order to discover the nature of the first 
equation in Equation (4-10), consider a small variation from a single van der Pol solution 


u(t): 


? 


2(t) = u(t) + €(z) oy 


where u(t) is a solution of the van der Pol equation, Equation (1-1). Substituting z into 
Equation (1-1) results in: 


i+€+e(1-(u+e) )(u+é)+u+e =0 (4-12) 
which when expanded and eliminating the higher order terms in é results in: 
fite(l—u?)ot+u}+é+e(1-w E+E =0 (4-13) 


where the term in braces is again the van der Pol equation and is identically zero, resulting 
in: 


E+e(1-w)E+E=0 (4-14). 
Equation (4-14) is the behavior of a small variation, €, away from the orbitally stable van 
der Pol limit cycle, therefore € must be orbitally stable. Comparing Equation (4-14) to the 
first of Equation (4-10), w and & have identical solutions, therefore w(t) must be orbitally 


stable. 


Understanding the stability of the first Equation (4-10) leaves only the behavior of the 
second variational equation: 


v—e(1-2B—u’)v+(1+42eA)v = 0 (4-15) 


to be characterized in order to understand the stability of Equations (4-1). Equation (4- 
15) will be referred to as the variational equation for the remainder of this thesis. 


SIGNIFICANCE OF THE ZERO MEAN DAMPING SURFACE 


Studying the stability of the variational equation requires characterizing the solutions of 
Equation (4-15) in the three dimensional parameter space { & , sA, €B }. Ultimately, it 
should be possible to find the surface in { € , sA, €B } space which corresponds to the 
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transition from stable to unstable oscillations of the coupled equations. Such an analysis is 
feasible, however it would be formidable to complete using analytical methods. Such an 
analysis using numerical techniques is underway by Reinhall and Storti[19,20]. 


In order to reduce the complexity of the problem, the parameter space may be reduced to 
two dimensions by fixing or specifying any of the three parameters. Such an analysis for 
the case of B=0, i.e. no damping between the coupled oscillators, is presented by Storti 
and Reinhall [26]. 


Another obvious two dimensional subspace for studying the variational equation is the 
Zero Mean Damping (ZMD) surface of the parameter space. The ZMD surface is the 
locus of values for the damping parameter ¢B, as a function of €, resulting in zero time 
average damping in the variational equation. The ZMD surface is a logical place to study 
the behavior Equation (4-15) because it marks the boundary between the region of 
negative mean damping, where all solutions to Equation (4-15) will grow over time due to 
the energy added by negative damping, from the region where damping is positive on 
average and stable solutions may be possible. By examining Equation (4-9), the parameter 
eA does not appear in the damping term, so the ZMD surface will be a function of € and 
€B only. The ZMD surface is found by requiring the time average of the coefficient of v 
to vanish in Equation (4-15): 


(1-2Ba —u?) = 0 (4-16) 


where the angled brackets indicate time averaging. Clearly Bzmp will be a function of 
since u(t), found in Chapter 1, is a power series in ¢. 


The balance of this thesis will be dedicated to characterizing the behavior of the variational 
equation, Equation (4-15), along the ZMD surface. 


CALCULATING THE ZMD SURFACE 
Performing the time averaging using the solution for u(t) found in Chapter 1 and is a 


tedious, but straight-forward process. Since Bzyp is not a function of time, Equation (4- 
16) can be rewritten as: 


1 1 
Ba. reer (4-17). 


Using the solution u(t) generated in Chapter 1, the ZMD surface can be calculated to the 
same order in € as u(t). The first few terms of Bzypp are: 
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Bap =-7- r+ + — - 
2 32 9216 10616832 152882380800 





(4-18). 


A Padé approximant to the ZMD surface through O(e”’)is given by: 


2 —0.5 ~ 0.185€ — 0,0669e4 - 0.0110€° — 0.00132 «8 - 0.0000589«!° 
Mw T+ 03072 + 0.1174 + 0.0163 + 0.002028 + 0.0000654€!0 (4.19) 


where the fractions which actually appear in the Padé approximant have been replaced 
with 3 decimals of precision. Figure (4-1) is a plot of the Padé approximant [10/10] of 
Bzmp calculated using exact coefficients. 





Bzmp 


Figure (4-1) Padé approximant to the 
Bzup Vs. c. 
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Figure (4-2) Padé approximant to the 
ZMD surface in 3-D parameter space. 


Figure (4-2) shows no additional information from Figure (4-1), but gives a physical feel 
for the parameter space. The stability transition curves obtained in Chapter 5 should be 
plotted along this surface in 3-space. 


FLOQUET THEORY AND STABILITY TRANSITIONS CURVES 


Equation (4-15), the variational equation, is a Hill’s equation, meaning it is a linear 
differential equation with periodic coefficients. This has two important consequences in 
relation to examining the stability of the solutions of Equations (4-1) along the ZMD 
surface. 


Since Equation (4-15) is a Hill’s equation whose periodic coefficients have period 7, all 
stable, periodic solutions will have period x or 27. Recalling that u(t) is periodic, with 


period 2x, u(t) must have period 7; then the claim follows directly from Floquet’s 
Theorem [16]. 
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The second result of the variational equation being a Hill’s equation, is that the transition 
from stable to unstable behavior in the parameter space will always be through a stable, 
periodic solution. In combination with the first result, that stable, periodic solution must 
be of period x or 2x. To prove this claim, consider Floquet’s Theorem[ 14]: 


The regular system x = P(t)x, where P is an n X n matrix function with 
minimal period T, has at least one non-trivial solution x = y(t) such that 
x(t +T) = p(t) where p is a constant. 


The variational equation can be written in the form required by Floquet’s Theorem by 
selecting: 


V 
K = tt (4-20) 


and identifying: 


0 l 
i oe £(1—- 2Bay — r ay 


which has period = due to the presence of u(t). Therefore Equation (4-15) meets the 
conditions of Floquet’s Theorem. Now consider the possible values for up. If w>1 the 
solution x(t) becomes unbounded over time since the solution is multiplied by for each 
interval of time 7; similarly if u<-1 the solution will also be unbounded over time. In both 
of these cases the solutions to the variational equation are unstable. If -1<u<1, the 
solution 1s stable and decays to the origin over time. That leaves two possible values for 


Ww: w=1 or p=-1. 


If p=1, the solution is periodic with period 2, the same period as P(t). If p=-1, the 
solution is also periodic but has period 2. Identifying where the characteristic number, , 
takes on the values of | and -1 marks the transition between stable and unstable solutions 
of the variational equation. 


The fact that the transition between stable and unstable behavior is always through a 
stable, periodic solution with period x or 27 focuses the study of this thesis. Locating the 
periodic solutions of Equation (4-15) on the ZMD surface identifies the stability transition 
curves, where the stability transition curves are the curves which divide stable and unstable 
behavior of Equations (4-1) in the parameter space along the ZMD surface. Hence, the 
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next task is to identify where periodic solutions (of period m or 2m) of Equation (4-15) 
exist on the ZMD surface. 





CHAPTER 5: PERIODIC SOLUTIONS OF THE VARIATIONAL EQUATION 


As discussed at the end of Chapter 4, identifying the periodic solutions of the variational 
equation along the Zero Mean Damping (ZMD) surface is equivalent to finding the 
stability transition curves for the coupled oscillators modeled by Equations (4-1). The 
next task is to find those periodic solutions. 


SOLUTION METHOD 


The periodic solutions of the variational equation will be identified using Lindstedt’s 
method [14], which was also used in Chapter 1 to solve the van der Pol equation. The 
variational equation along the zero mean damping surface is restated here for easy 
referencing: 


¥ —e(1-2B—u?)v+(1+2eA)v =0 (Cat) 


Some characteristics of the Equation (5-1) are worth mentioning. The most significant 
feature is that unlike the van der Pol equation itself, the variational equation is a linear 
ordinary differential equation. Since u(t) can be determined to a suitably large order in 
as Outlined in Chapter 1, and Bzyp is the expansion for the ZMD surface, all the 
coefficients in the variational equation are known except A. As a consequence of being a 
linear equation, solutions to Equation (5-1) will also obey the principle of superposition. 


Even though Equation (5-1) is a linear differential equation, Lindstedt’s method is used to 
find v(t) because of the presence of u(t) which is a power series. In keeping with 
Lindstedt’s method, the first step is to rescale time such that t = w(e) t. Substituting this 
rescaling into Equation (5-1) and changing the differentiation to the new time scale: 


v" —e(1-2Byp —u’ ov’ + @?(1+2eA)v = 0 (5-2) 


where primes represent differentiation with respect to the rescaled time t. Notice that here 
the function @(€) was already determined in solving Equation (1-1) and is given by 
Equation (1-19). The frequency expansion must be the same as for the original van der 
Pol solution because the solutions to Equation (5-1) must have the same period or twice 
the period of its coefficients; this was discussed at the end of Chapter 4. Plugging in 
power series expansions for @, v, Bzmp, u and €A according to: 
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eA 


i=1 
n,2 


Bryoi= Bo 1 ee 
u(t) = u, (t) + Due (5-3) 


v(t) = v,(t) + > v,(t)e’ 


and collecting powers of € gives a series of linear, ordinary differential equations, the first 
three of which are presented here: 


v,+(1+2A,)v, =9 
v’+(1+2A,)v, =—-2A,V, + Vv, -2B,vi -— us vi (5-4). 


v7 +(14+2A,)v, =-2A,V, —2A,v, —2u,u,Vv, + Vv; —2B,vi — ug Vv; —2@,uU7 
Notice that the left hand side of each of Equations (5-4) is of the same form, namely: 
LHS = v/+(1+2A,)v, (5-5). 


Insisting that the solution of Equation (5-1) be periodic with period 7m or 27, as discussed 
above, the coefficient for v; in Equation (5-5) must be an integer squared. This follows 
directly from the Floquet theory discussed at the end of Chapter 4; the solution to 


Equation (5-5) is Cos(t,/1— 2A, ) which can produce v(t) with period x or 27 only if 
{1—2A, is an integer; therefore: 


al 
14+2A,=k? or A,= —z— where Keno. (5-6). 


This gives rise to an infinite set of starting values for Ao, and possible stability transition 
curves out of each starting value. The first few Ao’s are given here: 


A ee Ce a (5-7) 
0” De ye ana 9° eye 
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DEALING WITH Ao=-1/2 


When Ao = -1/2, the variational equation changes form significantly; the zeroeth order 
equation of Equations (5-4) becomes: 


Vio) (5-8) 
which has solution: 
v, — ©, +C,,t (5-9) 


However, since the goal is to find periodic, stable solutions to the variational equation, Co, 
= 0 so Vo = Co 1s the solution. Plugging in this vo and examining the first order equation 
yields: 


v'=-2A,C, (5-10). 


In order for v; to be bounded, either Cp = 0 or Ay = 0. The former yields a trivial order 
zero solution, but the latter is a valid method of solving Equation (5-8) with a bounded 
solution. By similar arguments used about Equation (5-9) it can be seen that v; = C). 
Following this pattern, the same equation arises at each order of €, so that the stability 
transition curve from Ao=-1/2 is in fact just eA= -1/2. 


SOLVING FOR THE STABILITY TRANSITION CURVES 


For the remainder of the possible values for Ao listed in Equation (5-7), the solution 
procedure is a straight-forward application of Lindstedt’s method. Whereas the solution 
of the van der Pol equation in Chapter 1 resulted in fixing the frequency of the response as 
a function of €, solving the variational equation for v(t) will give the coefficients for the 
eA expansion. 


Unfortunately, unlike the arguments made in Chapter 1 for the van der Pol equation, there 
are no obvious symmetry arguments to use in the solution of the variational equation. 
While the solution of the variational equation shows interesting patterns for some values 
of Ao, the overall pattern is not easily discernible. This results in a very costly 
computational procedure for determining the stability transition curves. 


The stability transition curves calculated below and Tabulated in Appendix E were 
generated using the Mathematica™ code found in Appendix B. One feature of the code is 
worth mentioning here. In an effort to make solving the second order differential 
equations more efficient, a specific solving routine called MySolve[ _ ] was written. This 
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routine takes advantage of the fact that at every order of ¢, the differential equation 
encountered always takes the form: 


v'+k?v, = f(Sin(mt),Cos(mt)) where m= 1, 2, 3,...,2i+1 (5-11) 


and allows the code to produce solutions without calling a general differential equation 
solver. The specialized solver significantly increases the computational efficiency. 


CONDITIONS FOR SUPPRESSING SECULARITY 


Following Lindstedt’s method and substituting the power series of Equation (5-3) into 
Equation (5-2) and collecting terms of the same order in ¢, the first few equations can be 
obtained: 


k? v(0, t] + v0, t] = 0 
k? v([1, t] +v(1, t] = —2a{2] v(0, t] +2 v[0, t] — 4 Cos(t]? v0, t] 
k? v[2, t] + v[2, t] = -2 a{2] v[0, t] —2 a{!] v[1, t] — 3 Cos[t] Sin(t] v0, t] + 


l 
Cos[t] Sin(3 ¢] v0, t] +2 V1, t] —4 Cos(t]? v"[1, 1] + = v0, 1] 
(5-12) 
where k* = 1 + 2 Ao as defined in Equation (5-6). The notation used in Equations (5-12) 
is consistent with the Mathematica™ code found in Appendix B. In this notation, v[i,t] is 


the same as v;(t), afi] is Aj, and the numbers inside parenthesis represent differentiation 
with respect to the variable t. 


Noting that the left hand side of each of Equations (5-12) are identical, it can be seen that 
the homogeneous solution for each order of € will be similar: 


vali, t] = ci] Cos(kt] + sfi] Sin{k t] (5-13) 


where c[i] and s[i] will be used consistently to represent the coefficients of the 
homogeneous solutions of the i" order equation. The particular solutions will not be as 
consistent, but will be of the form: 


2i+l 


Vpli, t] = » (cli, m] Cos{m kt] + s{i, m] Sin[mk t] } 
n-0 (5-14). 


The c[i,m] and s[i,m] are tabulated in Appendix D. 
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In the course of solving Equations (5-11), two unique solutions for suppressing secularity 
will arise for each Ay [20]. The outcome of these two possible solutions for v[t] and A is 
that for each starting value of Ao, other than -1/2, there will be two curves along the ZMD 
surface which have valid, periodic solutions, with period m or 2x. These curves intersect 
at eA = Ap» and resemble a _ horn with its point downward and will sometimes be referred 
to as the “stability horn” arising from each Ay. Extending this result into the three 
dimensional parameter space of Figure (4-2), the two dimensional surface where periodic 
solutions to Equation (5-2) exist will resemble cones will their points coming out of each 
Ao. The regions within the cones are unstable, that is the solutions to the variational 
equation will grow over time, making the origin of the variational space an unstable 
equilibrium point. The regions outside of the cones and above the ZMD surface will be 
asymptotically stable; all solutions will decay towards the origin of the variational space 
over time. The actual surfaces of the cones are configurations which yield orbitally stable 
solutions. 


One of the main difficulties in solving Equations (5-12) will be to identify this split in the 
stability conditions. The conditions are different for each Ao and will be considered one at 
a time. 


CASE 1; Ap = 0 


For the case of Ap = 0, Equation (5-6) yields k=1. Solving Equations (5-12) and 
suppressing the secular terms appearing on the right hand side of the first order term 
gives: 

—c{0] -—2a{1l] s[0] =0 

—2a{1]c[{0] —s[0] =0 (5-15) 


which has two possible solutions for a[1] and s[0]: 


{{sto] + -c{0}, aft] > —}, {s{0} + cf0], aft] > -=}} 
2 2 (5-16). 

The first solution will yield a curve to the right of the second when viewed on the ZMD 

surface. From this point on, the solutions for v[t] and ¢A are unique and the stability 

transition curves are given by: 


mecca sc (36 247 & 11657 «7 1224811 & 
EAgw=-st tatoo t + —_—_—_—_———— 
2 8 32 384 1024 442368 21233664 2548039680 
2 6 Teé 36 247 & 11657 €’ 1224811 2 
Ay = (See. ee EE Oe) 
2 8 32 384 «+1024 442368 21233664 2548039680 
(5-17). 
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The power series given by Equation (5-17) were calculated out to O(e*’) and converted 
into Padé approximants which are plotted tn Figure (5-1) with the power series. 






«——_——. Padé Approximants 





Power Series 


= i 2 3 4 


Figure (5-1) Power series and Padé approximants 
[16/16] to transition curves out of A =0. 


CASE 2: Ay = 3/2 


For Ao = 3/2, Equation (5-6) gives k=2. Suppressing the secular terms in the first order 
equation of Equation (5-12) yields: 


— 2 a{l] s[0] =0 
—2 a{l] c[0]=0 (5-18) 


which requires a[1]=0 for a non-trivial solution. Suppressing secularity in the second 
order equation: 


osu u eae} sol S10 
c[0] 
—- — -?2 al? 0; =0 
3 a[2] c[0] (5-19) 


which has three possible solutions: 


{{al2] ~ -—, s{0] 7 0}, {al21 > > c[0] > 0}, {s[0] > 0, c[0] > 0}] (5-20). 
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The last solution is the trivial solution, but the first two again provide the split for the 
lower and upper stability curves. From this point, the solutions are unique and the stability 
transition curves are given by: 


3 & 1094 34196 297305 & 
éEAiy = —- — +o + OC) 


2 6 £3456 1990656 573308928 
3 2 534 #69836 740213 é 

CA = Eg as se 0 «!°) 
2 3 3456 1990656 2866544640 Cae 


The Padé approximants for the power series of Equations (5-21) calculated to O(e*’) are 
plotted in Figure (5-2) with the power series. 
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Figure (5-2) Power series and Padé approximants 
[12/12] to transition curves out of Ag=3/2. 
CASE 3: Ap = 4 


For Ao = 4, Equation (5-6) gives k=3. Suppressing the secular terms in the first order 
equation of Equation (5-12) yields: 


—2 afl] s{0] =0 
~2all] c{0] =0 (5-22) 


which requires a{1]=0 for a non-trivial solution. Suppressing secularity in the second 
order equation: 





9 s[0] 
e —2a{2] s[0] =0 
9c¢[0] 

= a —2a{2] c{0] =0 
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(5-23) 
which requires a[2]=-9/32 for a non-trivial solution. Suppressing secularity in the third 
order equation: 


5 c[0] _ 
a 2 a{3] s[0}] = 0 


= Zeal cl0] +S = 


(5-24) 
which has two possible solutions for a[3] and s[0]: 


{{s0] > -c{0], a[3] > =} {s{0] + {0}, af3] —}} 


(5-25). 
These two solutions provide the split for the lower and upper stability curves. From this 
point, the solutions are unique and the stability transition curves are given by: 


7 _ 4-28 5e , 28774 
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(5-26). 
The Padé approximants for the power series of Equations (5-26) calculated to O(e"”) are 
plotted in Figure (5-3) with the power series. 
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Figure (5-3) Power series and Padé approximants 
[16/16] to transition curves out of Ap=4. 
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CHAPTER 6: SUMMARY AND FUTURE WORK 


In Chapter 5, the first three pairs of stability transition curves along the Zero Mean 
Damping (ZMD) surface were determined. Combining Figures (5-1,2 and 3) yields the 
stability regions for the smaller values of A. In Figure (6-1), U marks regions that are 
unstable and S marks regions where orbitally stable solutions to the variational equation, 
Equation (4-15), exist. The stability transition curves are approximated by Padé 
approximants. 


Figure (6-1) Stability regions for coupled 
VDP oscillators on the ZMD surface. 


The three dimensional version of the same graph is given in Figure (6-2) where the three 
dimensions are the three parameters of Equation (4-1), {e, eA, B}. This is a graph of the 
stability transition curves found in Chapter 5 plotted on the ZMD surface of Figure (4-2). 


SUMMARY 


Using the techniques outlined in Chapter 5 and employing a symbolic mathematics 
package (Mathematica™) the first three pairs of stability transition curves for a pair of 
weakly coupled van der Pol oscillators were determined. All calculations were carried out 
exactly, that is exact fractions were used in all computations with no numerical roundoff. 
Using this analytical approach, the six stability transition curves seen in Figures (6-2) and 
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(6-2) were calculated to O(e*’). The power series coefficients for the transition curves are 
tabulated in Appendix E. 


FUTURE WORK 
There are three main areas for future work on this project. 


The first is to extend the range of eA covered in this thesis, namely find the stability 
transition curves coming from the rest of the Ap values given in Equation (5-7). Work on 
this area, using the code written for this thesis, is already in progress. The main difficulty 
in finding the stability transition curves coming from the next few values of Ao is showing 
that the curves are unique. Present conjecture is that there are more conditions needed to 
fully specify the upper curves than used in Chapter 5. To find what those conditions might 
be, research is underway to use Hill’s determinants to find the periodic solutions rather 
than Lindstedt’s method. 


Another area for future work was hinted at in Chapter 4 when the ZMD surface was 
introduced. The full challenge for this type of stability analysis would be to find the 
stability transition surfaces in the three dimensional parameter space. Currently, numerical 
calculations of this project are under way by Storti and Reinhall (REFERENCE>>>). The 
logical extension of this project would be to find the transition surfaces analytically to 
avoid numerical roundoff errors. 


The final, and ultimate continuation of this work is to physically realize the coupled 
relaxation oscillators modeled by Equation (4-1). Once these oscillators are accurately 
predicted by the mathematical models, they could be used to generate control signals for 
ambulatory robots, mechanical fish or similar mechanisms. The goal of this thesis was to 
further understanding toward reaching this goal. 
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Figure (6-2) Three dimensional plot of the 
stability transition curves along the ZMD surface. 
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APPENDIX A: MATHEMATICA™ CODE FOR VDP EQUATION 


This first block of code produces the first three terms of the of the VDP equation. 
Comments are included to help clarify the procedure. A second set of code is included 
which will start with a certain number of terms already in the solution and calculate two 
more terms. That code can be used to extend the solution as far as the user wants or until 
the computer runs out of memory. 


STARTER PROGRAM 


SLine = ©; 
SetDirectory["c:\\vdp\\canplex\\cdata"] 


(x pmult{f,g,x,nj} multiplies power series f by power series g 
where x is the variable, and keeps only through order n «) 


nnriies. 
pmult(f£_,g_, x Symbol, n Integer] := eae Coefficient[f, x, i] Coefficient[g, x, 3] 
i=0 420 


(x getX[n} creates a power series in e of order n with argument x[i,t] «) 


n 
getX[n_] := Expand[>\e*x{i, tl] 
is0 


(* xeven{] and xodd[] assign x[i,t] to be a complex cosine or sign «) 
2i+l1 
xeven[i_] := » Al} (Z? + 2) 


+ 


jal 
A J=a2 


2i+1 
xodd(i_j:= >\ IAi, 31 (2-2?) 
jal 


A Jj=2 
(x getW[n} creates the frequency power series «) 


n * 
getW[n_] := 1+ » e w[i] 
in2 


Aisc2 
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(x noslope[] enforces the boundary condition in Equation (1-1) «) 


noslope[i_] := Module[ {temp}, temp = Coefficient|[X[t] , e, i]; 
Solve[ (myD[ temp] /. Z—> 1) == 0, Afi, 1]]] 


(x myD[] and myD2[] are specialized derivatives «) 
myD[f_] :=f/.{Z-+IaZ,Z+12} 
myD2[£_] := myD[myD[f] ] 


ClearAll[x, w, X, W, temp, lhs, rhs,A, B, 1]; (« start with clean slate «) 


n= 3; (x order of solution «) 
X[t] = getX[n] ; (x build power series ine «) 
W= getW[n] ; 


(x build the right hand sides of Eq 1-5 «) 

rhs = pmultje, pmilt(W, prmilt(doc X[t] , 1 - pmilt[X[t] , X[t],e,n-1],e,n-1], 
e,n-1],e,n]; 

lhs = prmult(pmult(W, W, e, nj -1, d:¢,2; X[t] , e, nj; 


(* assign x[{1i,t] fomm to match Eq 1-15 «) 
Table[x[i, t] = xeven[i], {i, 0, n, 2}]; 
Table[x[i, t] = xodd[ij, {i,1,n, 2}]; 


(* back to building right hand sides «) 

mid = rhs- lhs /. Table[x©” [i, t] + myD[x[i, t]], (i, 0, n}]; 
mid= mid/. Table[x®” [i, t) + myD2[x[i, t]], (i, 0, n}]; 
X(t) = X[t] ; 


(x pull out each right hand side for each order of e «) 
temp = Table[Coefficient[mid, e, i], {i1, 1, n}]; 


(x free up some memory «) 
ClearAll[{lhs, rhs, mid, x]; 
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(* Solving the 1st order equation «) 
1= aes 


right = Collect| Expand[| ——— 


sec = ee Z]J]; (* Identifies the secular temms «) 

solRule = Solve[sec == 0, A[O, 1]]; (* Suppresses secularity +) 

A[O, 1] = A[O, 1] /. solRuleg3] ; 

right = right /. solRule[3] ; 

Coefficient[right, Zz Jet 
1- (23+ 1)2 

A[i, 1] = A[i, 1] /. noslope[i] {1] ; (* Enforces boundary condition +*) 

tempi] = 0; 


PED | , teble[z?*, (3, 0, 43] 


Table[A[i, 23+1] = , (J, 1, 1)]; (* Enforces Kg 1-15 %) 


(x Solving the 2nd order equation «) 
i1=2; 

right = Collect[ Expand|tempfill] , Table[z*7**, (5, 0, i}]]; 

solRule = Solve[Expand[Coefficient[right, Z]] == 0, wliJ]; (* Id’s the secular temms «) 
w[i] = w[i] /. solRule[l1j; 

right = right /. solRule[1] ; 

Coefficient[ right, 275+*] 


1-124: D? , (3, 1,i}]; (* Enforces Eq 1-15 +) 
= + 


Table[A[i, 23+1] = 


temp[i] = 0; 


(x Solving the 3rd order equation «) 
3; 


right = Collect| Expand| ———— 


solRule = ah ci sae Z)] == 0, Ali-1,1]]; (* Id’s the secular tezms «) 

A[i-1, 1] =A[{i-1, 1] /. solRule[{]]J; 

right = right /. solRule[1]; 

Coefficient[right, 27/**] 
ate 1) 

Ali, 1] = A[i, 1] /. noslope[i] [1] ; (x Enforces boundary condition «) 

tempi] = 0; 


oe ], Table[z?2**,, (5,0, i}]]; 


Table[A[i, 23+1] = , (3, 1, i}]3 (* Enforces Eq 1-15 «) 
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(x Saving the data... x) 
DeleteFile["afile.txt"] 
Save["afile.txt", A] 
DeleteFile|"wfile.txt"] 
Save ["wfile.txt", w] 
DeleteFile|["nfile.txt"] 
Save["nfile.txt", nj] 


(x This section puts the coefficients back in a 

form to use with sines and cosines...«) 
Table[Table[B[i, j] = 2A[i, J], {j, 1, 21+1, 2}], {1, 0, n, 2}]; 
Table[Table[B[i, j] = -2A[i, 3], {j3,1, 2i+1, 2}], {1,1,n, 2}]; 
Save["bfile.txt", B] 


(x Rebuilding X[t] in terms of sines and cosines «) 
2isl 
Table[x[i, t] = >| Bli, 3] Cos[jt], {i, 0, n, 2}]; 
jel 
A j=2 
2is1 
Table[x[i, t]= >| Bli, 3] Sin[jt], (i, 1, n, 2}]; 
j=1 
A j=2 


X[t] = Sum[x[i, t] e, (i, 0, n}|; 
(x Check out the work... «) 


Ww 
X[t] 


CONTINUATION PROGRAM 


The continuation program uses the same functions defined above, and the code 1s identical 
down through defining noslope[ | then change the code as follows: 





ClearAll[x, w, X, W, temp, lhs, rhs, A, i, nj; 
n<< "nfile.txt'"; 

old=n; 

n= old+ 2; 

Print["Starting, old = ", old," new = ",n) 
t1 = TimeUsed|] ; 

X[t] = getX[n] ; 

W= gew/[n] ; 

rhs = 1- pmult[X[t], X[t], e, n-1]; 
temp = pmult[W, %X[t],e,n-1]; 

rhs = pmult[rhs, temp, e, n- 1); 
Clear[tenp] ; 

rhs = rhs/. Table[e* — a, {i, old, n- 1} | ; 
rhs = rhs/.e~- 0; 

rhs = rhs/.a-e; 

rhs = pmult[e, rhs, e, n] ; 

lhs = pmult(W, W, e, nj; 

lhs = lhs-1; 

lhs = pmult[lhs, o;¢,2) X(t], e, n); 

lhs = lhs/. Table[e'a’*, {i, old+1,n)}]; 
lhs = lhs/.e- 0; 

lhs = lhs/. a~ e; 

Table[x[i, t] = xeven[i], {1, 0, n, 2}]; 
Table[x[i, t] = xodd[i], {1, 1, n, 2}]; 


mid= rhs- lhs /. Table[x®” [i, t] + myD[x[i, t]], (1,0, n)]; 


mid= mid/. Table[x (i, t] + myD2[x[i, t]], (2, 0,n)]; 
Print["updating X[t]..."] 

X[t] = X[t] ; 

Print[ "building temp..."] 

tenp = {mid/. {e"** .1, e+ 0}, mids. {eM 1, e+ 0}}; 
ClearAll[lhs, rhs, mid, x] ; 

A«<< "afile.txt"; 

w << "wfile.txt"; 
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i= old+ 1; 
ths = Collect [Expand temp[1}], Table[z*7**, (3, 0, i]; 
solRule = Solve[ Expand[Coefficient[ rhs, Z]] == 0, w{i]]; 
w(i] = wi] /. solRulefly ; 
rhs = rhs/. solRulefl]; 
Coefficient[{rhs, 275+] 


Table[A[i, 23+1] = 124+ 12 


»{j, 1, i}|; 


temp ly = 0; 
= old+ 2: 


rhs = Collect | Expand —=— 


solRule = Solve[ Expand| Coefficient|rhs, Z]] == 0, Afi-1, 1}]; 

Aji-1,1] =A[i-1, 1] /. solRulef]l]; 

rhs = rhs/. solRulef]] ; 

Coefficient rhs, 275+*] 
1- (23+1)2 

Aji, 1] = Afi, 1] /. noslope[ij ql]; 

temp{2] = 0; 

DeleteFile[{"afile.txt", "wfile.txt", "nfile.txt"}] 

Save["afile.txt", Aj 

Save["wfile.txt", w] 

Save["nfile.txt", n] 


ee), Table[z*7**, (3, 0, i}]]; 


Table[Aji, 23+1] = (3,1, i)]; 


From here the last few routines in the starter program can again be used to translate the 
solution back into sines and cosines and save the results in that form. 





APPENDIX B: STABILITY TRANSITION CODE 


The first program 1s “Recover.nb” which converts the coefficients of the complex VDP 
solution into real coefficients for use in the stability routines. It also calculates the ZMD 
surface coefficients. It updates these coefficients in the c:\stabil\data folder as well as the 
frequency coefficients. This code needs to be run only once after generating as many 
terms in the VDP solution as necessary. 


RECOVER.NB 


Directory structure 1s assumed to be 


C. 
-complex THIS HAS THE VDP SOLUTION STORED IN IT 
- stabil DATA FOR THE STABILITY CURVES 
SLine =O, 


SetDirectory['"'c:\\vebp"'7] ; 
n<< "camplex\\cdata\\nfile. txt"; (x mutber of terms available in VDP solution «) 


Reads in solution to complex solution and converts them to coefficients for real solution 
A<< "camnplex\\cdata\\afile. txt"; 
Table[Table[B[i, 3] = 2A[i, 3], (3,1, 2i+1, 2}], {i, 0,n, 2}]; 
Table[Table[B[i, Jj] = -2A[i, Jj], (3,1, 2i+1, 2}], {i,1,n, 2}]; 
DeleteFile|["canplex\\cdata\\bfile.txt"}] ; 
DeleteFile; "stabil\\data\\bfile.txt"] ; 
Save["canplex\\cdata\\bfile.txt", B] ; 
Save["stabil\\data\\bfile.txt", B] ; 
The next four lines calculate the coefficients of the ZMD surface 
2i+k 
Table[x[il = >' = Bi, gi (2° + Zi (1, Onn, 2)]; 
jz 
A j=2 
2i+1 1 : 
Table[x[i] = }' “> IBLi, 3 (29-25), (i, 1,n, 23]; 
j=1 
A je2 
A 1 : : 
Table[b[i} = DL ( Expand] -— x[3] x[i-3]] /-2-- 0) pi, OF De) 
4=0 
b[ 0] = b[O] + =; 


Saves the ZMD coefficients, updates the frequency coefficients in the stabil folder and 
updates how many terms are available in the VDP for the stability routines 





SO 
DeleteFile|'"stabil\\data\\crit.txt"] 
Save["stabil\\data\\crit.txt", b] 
DeleteFile|["stabil\\data\\wfile.txt"}] 
CopyFile|["canmplex\\cdata\\wfile.txt", "stabil\\data\\wfile.txt"] 
DeleteFile["stabil\\data\\nfile. txt") 
Save["stabil\\data\\nfile.txt", n] 


The next program will start the solution for the stability transition curve: 


STARTER.NB 


Assumes the same directory structure as Recover.nb 
SLine = © 
SetDirectory| "c:\\vdp\\stabil\\data"] 


pmult[ ] is the same as in the VDP code. getV[ ] 1s equivalent to getX[ | from VDP and 
getA[ ] makes the power series for the transition curve. 
n mi 
pmlt(f_,g_, x Symbol, n Integer] :- }' )'x'*) Coefficient(f, x, i] Coefficient(g, x, 3] 
isO 320 


order 
getV[order ] := >, e"v[n, t] 
n=0 


onder 
getA[order ] := » e"a[n] 
n=O 


update[ | updates the RHSs of Equation (5-12) as each v[i,t] is determined. 
update[temp_, i_] := Module[{ttt), ttt = tenp; 
tttfi- 1} =-0; 
Table[ttt[iJ = ttt{i] /. vij, tl > vij, tl, (3,0,i-1)]; 
Table[ttt{il = ttthiy /.v°' (5, t] > devi, tl, (3, 0, 1-13]; 
Table[ttt{il = tttgig /. v5, t] + o¢,2) VL, t], (5,0, i-D]; 
ttt] 


Specialized ODE solver discussed in Chapter 5. 
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mysolve[i_, right _] := Module {n, temp} , 
n=2 (1+ 1) + 1: 
temp = C[1] Cos[k t] + s[1] Sin[k t] ; 


X1 Ccoefficient[right, Cos[mt] ] Cos[mt] 
temp = temp+ i a a_i ; 
“, Coefficient[right, Cos[mt]] Cos[mt] | 
temp = temp+ a a ; 
Rel Coefficient[right, Sin[mt]] Sin[m t] 
temp = temp+ a =77 ce, an ; 
- . y' Coefficient[right, Sin[mt]] Sin[mt] | 


-m2 2 


temp] 


Build the power series for the frequency, W, the VDP solution, U, and the ZMD surface, 
B. 
= af 
W-= 1+ e w[i] 7 
Ais2 
> 
U= > 2 uli, t]; 


i=0 


5 
B= e'b[i]; 
23 


A in2 
ClearAllin, v, V, A, a, temp, right, s, c] 
n= 5; (* how many terms to calculate «) 
V[t] = getv[n] ; 
A= getA[n] ; 
a[0O] = 0; (x which value of a[0] to go after «) 
k= -V1+2a[0) (» k of Eq (5-6) «) 


ths is the right hand side of Equation (5-12) prior to separating the powers of e. 

Once they are separated, each power is stored in temp[i]. 

rhs = pmult[pmult[ew, 0 V[t], e, n], 1- pmlt[U, U, e, n-1] -28, e, nj - 
2pmlt(A- a0, V[t], e, n) -pmilt(pmilt[W, W, e, n) - 1, Ot,2; Vt], e, n); 

temp = Table[Coefficient[rhs, e, i], {1, 1, n}]; 

Clear| rhs] ; 


Read in the coefficients for VDP, frequency and ZMD surface (generated by recover.nb). 
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u<< "newu. txt"; 
w << "wfile.txt"; 
b<< ‘Ere -txc'- 


Solve the zeroeth order equation, note the rhs is zero for this equation. Last line is the 
condition of unit displacement for the perturbation. It can be replaced with other 
conditions or left out entirely. 

1 ="0; 

v[i, t] = mysolve[i, 0]; 

v[i, t] = v[i, t] /. Solve[ (getv[i] /. t+ 0) == 1, cfi]] {1; 


Solve the first order equation. First update the rhs to reflect v[0,t], next find the secular 
term coefficients, s1 and cl. Next solve to eliminate the secular terms, in this case for a[1] 
and s[0]. Next update the coefficients and the rhs with the solution just found and finally 
call the specialized ODE solver to get v[1,t]. Again the last line imposes unit displacement 
and could be left out or replaced by other conditions. 

Note: For this case, a[O]=0, picking sol{f{1]] gives the lower curve while sol{{2]] would 
give the upper curve. 

i= be 

temp = update[ temp, i] ; 

right = Collect|TrigReduce|[tempgij], (Sin[t] , Cos[t] }]> 

81 = Coefficient[right, Sin[t]]; 

cl = Coefficient[right, Cos[t] ] ; 

sol = Solve[{sl1== 0, cl == 0}, {a[i], s{i- 1]}] 

afi] = afi] /. solgl]; 

v(i-1, t] =v{i-1, t] /. soll]; 

s[i-1] =s[{i-1]/. solglj; 

right = right /. soli; 

v[i, t] = mysolve[i, right] ; 

v{i, t] = v[i, t] /. Solve, (getv[i] /. t> 0) == 1, cfij] qld; 


Solve second order equation, just like above. 

i= 2: 

temp = update/[tenp, i] ; 

right = Collect[TrigReduce[tempgiy], {Sin[t] , Cos[t]}]; 
Sl = Coefficient[right, Sin[t]] ; 

cl = Coefficient/[right, Cos[t]] ; 

sol = Solve[{sl== 0, cl== 0}, {afi], s{i-1]}] 

afi] = a[i] /. solf{1]; 

s{i-1] =s{i-1] /. solf{1]; 

right = right /. solf{1}; 
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v[i, t]) = mysolve[i, right] ; 
v[i, t) =v{i, t] /. Solve{ (getv[i] /. t- 0) == 1, c[1])] ql]; 


Same thing for third order... 

i= 5 

temp = update[tenmp, 1] ; 

right = Collect[TrigRediuce[tempfi]], {Sin[t] , Cos[t]}]; 
sl = Coefficient[right, Sin[t]] ; 

cl = Coefficient[ right, Cos[t] ] ; 

sol = Solve[{{sl == 0, cl== 0}, {a[1i], s{i-1]}] 

a[i) = a[i] /. solgl]; 

s[i-1] =s[i-1] /. solf{l]; 

right = right /. soljl]; 

v[i, t] = mysolve[i, right] ; 

vii, t) = v[i, t] /. Solve[ (getv[i] /. t- 0) == 1, c[1])] 01]; 


Same thing for fourth order... 

1=4; 

temp = update[temp, 1]; 

right = Collect| TrigReduce|tempyij], (Sin[t] , Cos[t]}]; 
sl = Coefficient[right, Sin[t] ] ; 

cl = Coefficient{ right, Cos[t]]; 

sol = Solve[{sl == 0, cl == 0}, {a[i], s{i-1]}] 

a[i]) = a[i] /. solf{ 1]; 

s[i- 1] =s[i-1] /. solf{1]; 

right = right /. solf{l]; 

v{i, t] = mysolve[i, right] ; 

v{i, t]) = v[i, t] /. Solve[(getV[i] /. t~> 0) == 1, c{iJ] [1]; 


Same thing for fifth order... 

n= 5° 

temp = update[temp, i] ; 

right = Collect[TrigReduce|[tempyig], (Sin[t] , Cos[t]}]; 
s1 = Coefficient[right, Sin[t]]; 

cl = Coefficient[right, Cos[t]]; 

sol = Solve[{sl1 == 0, cl==0}, {afi], s[i-1]}] 

a[i] = a[i] /. solql]; 

s[i-1]) =s[i-1] /. solgl]; 

right = right/. solgqlj; 

v[i, t] = mysolve{i, right] ; 

v[i, t]) = v[i, t] /. Solve[ (getV[i] /. to 0) == 1, c[1]] [1]; 
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Save the results. a0l.txt 1s the lower curve out of a[OJ=0. nOl.txt is how many terms of 
v[i,t] have been solved. 
DeleteFile|['"a01l.txt"} 


Save[''a01l. txt", A] 
DeleteFile["v0O1l.txt"] 
Save["v0O1.txt", V] 
DeleteFile;/'"n01.txt''} 
Save[''n01. txt", nj 


CONTINUE.NB 


Continue will take the results of either starter.nb, or itself and extend the results by four 
terms. The functions down through and including mysolve[ | from starter.nb need to be 
included as well as MyExpand[ ] which uses Euler’s identity to update the right hand side 
of Equation (5-12) more quickly: 


MyExpand[mess_] := Module[{Z, temp}, 


temp = Expand| mess /. {Cos[m_t] > = (2"+ 2°"), Sin[m_t] > 5 I(z™-Z™), 
1 1 iL en 
Cos[t] - > (Z+ = , Sin[t] > "> = (2-—)}] : 


Expand|temp/. {Z™-~ Cos[mt] + I1Sin[mt], 2° + Cos[2t] + 1Sin[2t], 
Z—- Cos[t} + ISin[t}}]] 
The next section of code checks to ensure enough terms are available in the VDP solution 
and then builds the nght hand sides of Equation (5-12) much as above. It prints out lines 
to let the user know what order terms it is trying to calculate: 


n<< "nfile.txt"; 

m=n; 

Print[m, " terms in VDP solution available."] 

n<< "nO1l. txt"; 

min = n+ ie 

N= min+ 3; 

Print["Starting with ", min, ", going to ",n, "."] 
If[n>m, Print["Not enough terms! !"] ; Quit[], Clear[m}]] 





5S) 
n 
W= 1+ om e w[il]; 


i=2 
A i=2 


n 
U= euli, t]; 
n > 
ZMD = >a e'b[i]; 


1=0 
Ais2 


ClearAll[v, V, A, temp, rhs, s, c] 
V[t] = getV[n] ; 
A= getA[n] ; 
rhs = pmult[ew, 6, V[t],e, n] ; 
temp = 1- pmult[U, U, e, n] -2 2D; 
rhs = pmult[rhs, temp, e, n, min] ; 
temp = 2pmult[A- a[0], V[t], e, n, min] ; 
temp2 = pmult(W, W,e, nj -1; 
temp2 = pmult[temp2, 6;t,2) V[t] , e, n, min] ; 
rhs = rhs - temp - tenp2; 
ClearAll|[temp, temp2] ; 
temp = Table[Coefficient[rhs, e, i], {i, min, n}]; 
Clear| rhs] ; 
<< "a0l.txt"; 
<< "VOl.txt"; 
k= V1+2a[0] ; 
oie 
Table[uli, t] = >' Bi, j] Cos[jtl, (i, 0, n, 2}]; 
ke 
2 dal 
Table[uli, t] = >\ Bli, 3] Sin[jtl, (i, 1, n, 2)]; 


5-1 
A jJ=2 


B<< "bfile.txt"; 
w << "wfile.txt"; 
iec< “orit.txt"> 


The next section calculates the four new terms, and are just like the equivalent sections of 
starter.nb except for using MyExpand[ ] instead of Expand{[ ]: 


1= min; 

Print| "Working on term: ", i] 
temp = update[temp, i]; 

rhs = MyExpand|tempgi - min+ 1]] ; 





sl = Coefficient[rhs, Sin[k t]] ; 

cl = Coefficient[ rhs, Cos[k t] ] ; 

sol = Solve[{sl == 0, cl == 0}, {a[i] , s[1 -1]}]; 

a[1] = a[i1] /. soll]; 

s[i-1] =s[i-1] /. sol[l]; 

rhs = rhs /. soll]; 

v[i, t] = mysolve[i1, rhs] ; 

v[i, t] = v[1, t] /. Solve[ (getv[i] /. t> 0) == 1, c[1]] [1]; 


1 = min+1; 

Print| "Working on term: ", 1] 

temp = update|[temp, i] ; 

ths = MyExpand|tempyi- min+ 19] ; 

sl = Coefficient[rhs, Sin[k t]] ; 

cl = Coefficient|[ rhs, Cos[k t]] ; 

sol = Solve[{sl == 0, cl == 0}, {a[i], s[i-1]}]; 
a[i] = a[i] /. sol[l]; 

s[i-1] =s[i-1] /. sol[1]; 

rhs = rhs /. solfi] ; 

v[i, t] = mysolve[i, rhs] ; 

v[i, t] = v[i, t] /. Solve[ (getV[i] /. t- 0) == 1, c[i]] [1]; 


1= min+ 2; 

Print[ "Working on tem: ", i] 

temp = update[tenmp, 1] ; 

rhs = MyExpand|[temp[i- min+ 19] ; 

sl = Coefficient[rhs, Sin[k t]]; 

cl = Coefficient|[ rhs, Cos[k t]]; 

sol = Solve[{sl == 0, cl == 0}, {a[i] , s[i-1]}]; 
a[i] = a[1] /. sol[l]; 

s[i-1] =s[i-1] /. soljl]; 

rhs = rhs/. solji]; 

v[i, t] = mysolve[i, rhs] ; 

v[i, t] = v[i, t] /. Solve[ (getV[i] /. t- 0) == 1, c[1]] [1]; 


1= min+ 3; 

Print[ "Working on tem: ", i] 

temp = update|[temp, 1] ; 

rhs = MyExpand|tempji- min+ 19] ; 

sl = Coefficient[ rhs, Sin[k t]] ; 

cl = Coefficient[rhs, Cos[k t]]; 

sol = Solve[{sl == 0, cl == 0}, {a[i], s[i-1]}]; 
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a[i] = a[1] /. soll]; 

s[i-1] =s[i-1] /. soljl]; 

rhs = rhs/. solfl]; 

v[i, t] = mysolve[i, rhs] ; 

v[i, t] = v[i, t] /. Solve[ (getV[i] /. t- 0) == 1, c[1]] [1]; 


Finally save the results as above: 


A= getA[n] ; 

V[t] = getv[n] ; 
DeleteFile, "a0l.txt''] 
Save["a0l.txt", a] 
DeleteFile|[ "v01l.txt"] 
Save["v01. txt", v] 
DeleteFile|'"n01l.txt"] 
Save["nOl. txt", nj 
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APPENDIX C: SOLUTION TO VDP EQUATION 


The analytical solution to Equation (1-1) through order &’’: 
X(t) = )23@es +6 (ace + ae Cos[3t] - 28 Cos[5 t}] + 
16 96 


7 13 Gos (re) 7 47 Cos[3t] HOSS COS |S] - 2149 Cos[7t] 61 Cos[9t] 


—————— Ee ee mice 
12288 1536 “ 21648 110592 20480 
ef ( 6479 Cos[t] c 48437 Cos[3t] : 259945 Cos[5t] ¥ 


6635520 35389440 31850496 
A253767 Gos|7 t] 480523 Cos[9t] 4937537 Cos[ile) 715247 ¢os) 136) 


398131200 73728000 2654208000 3715891200 
| 377080601 Cos[t] | 36921629Cos[3t] 4094345 Cos[5 t} 


eee eed 
1712282664960 71345111040 9172942848 
107409503789 Cos[7t] i 188691979247 Cos[9t] 2067384798557 11 Cos 1G) P 
45864714240000 59454259200000 936404582400000 
9257 312351) Cosi is te) 7184260289 Cos[15t] 392636471 Cos[17 t] 
ee 
1310966415360000 4661213921280 29964 946636800 
el0 (_ 44898976356847 Cos[t] _ 4475878408049 Cos[3 t] 
30204666209894 40000 50341110349824000 
611229379671 Cos[5T) : 21294508407929 Cos[7 t] 7 1838562603094127 Cos[9t] - 
2958824445050880 165112971264000000 2621932830720000000 
1299787042060760957 Cos[11 t] 7 11188471201628174701 Cos[13 t] : 
1321454146682880000000 14800286442848256000000 
217736422934204291 Cos[15t}] 14852537 752791657 7 Gos) 17 
rns we 
78934861028524032000 1522315176978677760000 
1600597132693073 Cos[19 t] 29654422883 Cos[21 t] 


108736798355619840000 32288758824 960000 
E fee spi Sin[3t] +E es = Sin[3t] - is Sin[5t] + ul 
4 4 256 256 576 576 
5 is 12971 Sanalcte) 259 lesamy | 3 ts] 
4423680 294912 
52885Sin{St] _ 110621Sin(7t) | 7457Sin{9t) _ ee ; 
2654208 6635520 1228800 7372800 
Be | SS 114653 Sini tc) 7 82937 Sin[3t] 7 1939625 Sim? 5 t] F 1061235889 Sin[ 7 t] 
44590694400 212336640 764411904 191102976000 
179467921 Sin[9t] , 229342933 Sim| lit) a 195050323 Simjli3 c] ‘ 138697 Sin[15 t} 
35389440000 92897280000 346816512000 2774532096 
<2 | 9822191352131 Sin[t] f 89866136849 Sin[3t] 1702014109 Sin[5t] 





Sin[7t]] m 


251705551749120000 342456532992000 12328435187712 
2088126932941 Sin[7t] 41070940258801 Sin[9t] 1694102483205143 Sin[11t] r 


ye fea 
2751882854400000 24970788864000000 10487731 32288000000 
662514241768639 Sin[13t] 507370714969Sin[15t]  15080405867887 Sin[17t] 


ee 


734141192601600000 1740186530611200 302046662098944000 
eee 


13484225986€5600 
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Frequency Expansion through order 20: 
e 17Je! 35 &® 678899 e® 
W = Le met meen tee Ces + 
16 3072 4884736 5096079360 
28160413 e1° 16729607288111 e” 5722795344767278507 e!4 
a I oe 0 a app gp 
2293235712000 3698530556313600000  5219366321069752320000000 
1846779765852173498887007 e!° 320250638 6889338 212533493973243 el 
cs cf eee 
14731139504587268947968000000000 41577168137747107878744883200000000000 
4038372199485064617691118289043994731 e29 


3872464178 6152554301399954253 41 440000000000000 


NUMERICAL RESULTS 


These lists of numbers represent the coefficients of the cosine or sines of the VDP 


solution. c[2] for example shows that x2(t) = -0.125 Cos(t) + 0.188 Cos(t) -0.0521 
Cos(t). 


Cosine Coefficients through order 25: 


c[0] 

{Z.} 

c[2] 

(20: 125 e019" -0.0521} 

c[4] 

(0.00594, -0.0306, 0.0392, -0.0194, 0.00298} 

c[6] 

(0.000976, 0.00137, -0.00816, 0.0107, -0.00652, 0.00186, - 0.000192} 

c[8] 

{-0.00022, 0.000518, 0.000446, -0.00234, 0.00317, -0.00221, 0.00085, -0.000168, 
0.0000131} 

c[10] 

{-0.0000149, -0.0000889, 0.000173, 0.000129, -0.000701, 0.000984, -0.000756, 
0.000351, -0.0000976, 0.0000147, -9.18x10 7} 

c[12] 

(0.000011, -0.0000129, -0.000034, 0.000059, 0.0000364, -0.000214, 0.000313, 
-0.00026, 0.000139, -0.0000479, 0.0000104, -1.26x 10°, 6.56x 10°} 

c[14] 

{-3.84x 107’, 5.51«x10°°, -4.1x107®, -0.0000121, 0.0000203, 9.83 10°°, -0.0000661, 
0.000101, -0.00009, 0.0000533, -0.0000216, 5.91x10°°, -1.04x 107°, 1.06x 10°’, 
Agel O°} 

c[16] 





60 

{-5.56x10°8 1.1x10#, 2.13x10°°, =1.4x 106% =40o4emie: “7 0p lee 
2.46x 10°°, -0.0000206, 0.000033, -0.0000312, 0.0000201, -9.22x 10°, 3.x10°, 
5G. 79x 10°’, 1.01x10° , —-6e< 10. pceaemie 

c[ 18} 

(9.17x10-°,. -3:22«10°', =8. 11k 10:7, 7:855e10n ade 210; ieee 
2.43x10°°, 5.36« 10°’, -6.42« 107°, 0.0000108, -0.0000108, 7.52x 107°, 
gee8ix 10°, 1.42x10°°, -3.85 10°’, 7.42<10°°, -9.55% 10°, 11S50 10, 1 2 oo gee 

c[20} 

(PEAG 10) O84 x 10 (pn 22 10 | p= 1. 12x10) pe ee els > — 160 ome one 
Bree l0: esr x10", =2.Olel0 °, 3.5710", ssegoxde 712. 7G sch: eee te 
Gio 10°', =2.2 10°, 4.65% 10 °, -7.78% 100, 6282x110, 650510 7 ne or Oe 

c[22] 

peece 10) malGoc10.°, 1,41x10 °, =4549% 10°", 6.79% 10) 1. 02 atOn = 5202 10, 
BiG a > O4emlGn = 7.0410, -6.24n10 Ina 10 - Ieeralioae 
102 10; NOS 10) 2.76x 10), 9-72 10), 2.63x10 , =5es5 0, 

Wm 10°, =8.01x10~, 4.95x10 “, -1.49x10 7} 

c[24] 

fees. 108 3268.10, 6.03~ 10, 504x100, @1.6o~ 10 “pea.zeur; 
pac 10h a= Osan 5.590010» 10210 2 aaa tO) Soe 101, 

Breer 40a) 1, 3-74x10) , =2.36x10 |, 1 lexl0 , —da52~40) 
ee 10° =3.99~™x Wr, 5-92x10 , -7.84%10 , 7.17% 10 -, -4e03- 00°, 1-05. Gnas 


Sine Coefficients through order 25: 

s[1] 

(0.75, -0.25} 

s[3] 

{-0.0273, 0.082, -0.0608, 0.0122} 

s[5] 

{-0.00293, -0.00879, 0.0199, -0.0167, 0.00607, -0.00075} 

s[7] 

(0.000743, -0.000391, -0.00254, 0.00555, -0.00507, 0.0024, -0.000562, 0.00005} 

s[9] 

(0.000039, 0.000262, -0.000138, -0.000759, 0.00164, -0.00162, 0.000902, 
-0.000292, 0.0000499, -3.46x 10°°} 

s[11] 

{-0.000035, -9.57x 10°, 0.0000953, -0.0000478, -0.000227, 0.000505, -0.000527, 
0.000332, -0.000132, 0.000032, -4.32x10°, 2.45x10°7} 

s[13] 

{1.51x 10°, -0.0000106, -4.09x 107°, 0.0000335, -0.0000174, -0.0000683, 0.000158, 
-0.000175, 0.000121, -0.0000556, 0.000017, -3.3x 10°, 3.67x 10’, -1.76x 10°} 
s[15] 





6] 

{1.72x 107% 1.56% 107°, =4.06x 1075ys21-.73 1079 OeOGOdmieeeG >. oe 

-0.0000206, 0.0000501, -0.0000583, 0.0000437, -0.0000225, 8.14x10°, 

22202 10° 93. 26x 10 723.0810" pullmce lone 
s[17] 

(=2.6G.0105', 4.13107’, 6.46.10; = 11 Miex 10po onda ie ee slide 

2.4410, =6.16x 10°, 0.000016, -0.0000196, 0. 0000157, 8.9110, 

2G GvulOno eos 10° pee 10 |». - Sele 10 2 be len On imine 
s[19] 

eg 03 We i410’, 1.57% 10’, 2.5910 |, <5. 1610", — 2.2400 ine elem 
Bor lOmin eo 107°, 5.15 107°, ~6.62 10°, 5362. 10 °, = 3. 45.101, 

1S GscOre = 50c7s 10° wees 10"', 2. 4clge) Deol. 10 => 10, ee Oa 
s[(21] 

Pei Ie. 02410, = 554% 10°, 5-4 5am °, See 10. eeacu On 

57 AG 0e e459 Women = 3.4310" » =5.37« 10 ies x 10 eee ns 

DiGi. Or 13 see 10°, 6. 59a eu. 51 ese) a= ISOs 

2 Ax 10m 2.66 0mm 7S< 107, - 5.1510 >} 
s[23] 

ONO? < dlQmmno 52 x 105s, 45x 105 2.18 x 108, 1. 86x 103. 59x 10am 

EG SicalUmee— 2 4 om eel). 56% 10m 1.29 x 10m 1. 54x 108 es, 4 

=7.57 x lene 7. 15 Wome 5. x 107 2-69 107 el 13x 10°, 3eegs 10 °, 

PCnaAeelOm 119 10e 5x 10ers 2.4x10 ede ail 10) elon 
s[25] 

fe 1. 66mm - 4512108) 3.9110) 28.87 x 10, Seo 10 sone 0 eo), 
ED Daina — 7. 3yeln 4, 5.29% 0p » =4. 73 x 10 4. 31x10 plea 
PSG TU 2.54 lime o7x 10", 1.08x100, 4.9% 10° eon 

Bere cil0: ale 108 we 11.94% 107°; Dea 210) 2 ad 10, Teo 101 el 





APPENDIX D : SOLUTIONS TO THE VARIATIONAL EQUATION 


This appendix includes the solutions to the variational equation through ¢ for the cases 


covered in Chapter 5. Note that these solutions can be used to check if the code in 
Appendix C is working correctly. 


Curves out of Ao=0: 








a Cos[t] 1 3Sm[t] 1. 
Lower = Cos[t] + Sin[t] + c{- ; + a Cos[3 t] + Gia Sin[3 u) + 








»( 19Cosft] 1 5 7Sinft) 1 5 5 Cos[t] 
€ {-———+ — Cos[3 t] - —— Cos[5t] - + — Sin (3t) - < siniSu} + « i 
i228 192 1922 16 516 
7, cos _ TCost74} _ SSSinit] 53 89Sini5t) st) 
256 768 1152 2304. 768 2304 1152 
7 4249Coslt] _S7Cosi3t] 3371 CosfSt] _ 2615Cost7t) 61 Cos{9t]_ 439Sinft) 
552960 2048 110592 221184 40960 552960 
107Sin[3t] 1631Sin[St] 1885Sin(7t] 61Sim[9t]) _¢ / 1013 Cosft] 
9216. 110502 + +221184 ~~ ~—40960 J+< 7372800 
6707 Cos[3t)  40205Cos[St} 180211 Cos[7t]  40709Cos[9t] 5533 Cos[11t] 
211840 5308416. ~~«-26542080.~«~S=«*dA74S600..~—~»~=~«*2« 745600 
8311 Sinft]  21127Sinf3t] 91811 SinfSt] 319577 Sinf7t) _158179Sinf9t] $533 Sin[ 114 
22118400  —-2211840 5308416 -—-- 26542080 44736800 -—«14745600 








' Cos[t] | 3 Smn[t] l 
Upper = Cos[t] — Smt] + { 3 


- 8 Cos[3 t] a = 8 Sm([3 t] ate 





| 19 Cos[t] Ef Bq 5 Be 7 Smt] l Sin(3t 5 ~— sin(St)} o(- 5 Cos[t] 
_ ———— + — Cos _—— ——— as eS 
: 192 8 i 10 a 4 








576 
7 19 _ 2Cos(7t} 55Smit}_ 53 g9Sinf5t] 7 Sin{7t] 
sap Cos[3t]- = Cos[5t] + ——— - Sait —— ee 
256 768 1152 2304.” 768 2304 1152 
4( 4249Cos[t] 57 Cos[3t] 3371 Cos[5t] 2615 Cos{7t] 
| 552960 2048 (2110592  ~—«-221184 
61Cos[9t]  439Sinft]  107Sinf3t] 1631 Sim(5t] 1885Sin[7t) 61 Sm[9t] 
4090. 552900 9216. +  +4110592 221184 © 40960 J 
5( 1013Cos[t]  6707Cos[3t)  40205Cos[St] 180211 Cos{7t]  40709Cos[9t] 5533 Cosf11 
- 7372800 (2211840 +«~~=«5308416.~=S*é<“i«t*«éSDOBO.-—-”~”—C«~O*SNKTASGOOOD~=—~S~*~S~«*ASGCO 
8311 Sin{t] 21127Sm[3t] 91811 Sin[St) 319577Sin7t] 158179 Smp9t] 5533 Sin 114 
22118400 2211840 5308416. + +«-26542080 +~~«44236800+~=~S~S*«*2: 745600 


Curves out of Ag=3/2: 
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71 a 7 
Lower = Cos[2t] + 3 —— Cos[2t] + —— Cos[4t] - — Costs} + 
576 144 192 


4{ 15443Cos[2t]  2809Cos[4t]  1427Cos[6t] = 1127 Cos[8t}] 119 Cos[10t] 
Tp | + — | + 
( 1658880 82944 36864 69120 55296 
bike bw: aj i3sm(2t) 4/7 © 239 Sm[6 t] Sm 
(= Sin{2t]- — = Sint] +€ [—— + —— Sn[4t] - ———_——_ + —— sin'8t)} + 
24 1728 576 4608 576 
5 (- 5407 Smf[2t] = 112483 Sm[4t] 285299 Snf6t]  65987Sm[8t]) 22187Sm[l0t] = 1007 Sm[12t] 
Se ee enn rennet emery fe ne RR 
: 13271040 9953280 13271040 4147200 4423680 1843200 


11 1 
DF — ¢/— — Cos[2t] + — Cos/[4t 
ppa { 48 os[Zt] + D os 1) 


5( 35Cos[2t] 181 Cos[4t] 239Cos[6t] S$ Cos[8t] 

i ( sgn) ASG 1 152 
218417 Cos[4t]  378463Cos[6t] 1117 Cos[8t] 

19906560 26542080 ~~—~—«+129600 

© suroty + 2 Se Aa saz Sint6ti} + 

2 1152 288 

7 46007 Si2t]  4385Sin{4t]  4883Sin(6t) 1127 Sin(8t) 
3317760 165888 221184 


5( 319213Cos(2t] 

}. (- 79626240 
22187Cos[10t] 1007 Cos[12t] 
~ 8847360 3686400 J 


119 Sin[ 104] | 
138240. 110592 


Curves out of Ag= 


Lower = Cos[3t] — Sm[3t] + 











Ceres 3 3Sin(t} 359 3 af 329Costt 
c{ — — Cos[3t] — — Cos[5t] + + —— Snn[3t] - — — Sint 5t}} + (- rs 
8 16 16 240 640 

233Cos[3t] 479Cos[5t] 27 51 Sint] 73913 on ion 7 

ee a = Coe ~~ Sint} + — Sin 70 + 

1280 1280 640 128 57600 256 640 
,( 89813Cosft] 76151 Cos[3t] 139 67 Cos{7t] 
: (= 2 ee EE G5 ee 

153600 230400 600 2048 

469.Cos{9t]  831Sin(t] _ 82393Sini3t]  679Sw/St] S311Sin(71] 469 Sinl91 ) 
46080 2048 57000 -~S—=«<CSdHS*é“‘é‘«CSC‘SC*«*:‘C‘CKOVNSC*”?Y#CS 


4 ( 385217 Cos[t] 105987247 Cos[3 t] 
| eS a 


614400 ~-—-:309657600 
389277 Cos[5t]  30401Cos{7t] 327931 Cos[9t] 


547Cos[11t] 273059 Smt] 
+ i 
1638400 409600 11059200 215040 614400 
37616912581 Sm[3t]  18809Sm[St] 249571 Sm[7t] 27691 Sm[9 y 547 Sm[11 t] 
qx. ———~- + Orr > 
23224320000 327680 6144000 2211840 215040 
5 ( 42447078481 Cos[t] 33657149027 Cos[3t] 9222921259 Cos[5 t] 1330129 Cos[7 t] 
61931520000 92897280000 / 30965760000 39321600 
46069577 Cos[9t] 88301 Cos[11t] 89371 Cos[13t] 400303339 Sm[t] 82024682603 Sm{3 t] 
+ Oooo 5 Oe ere 
9289728000 20643840 137625600 825753600 46448640000 
13130773 Sm[S5t] 16057111 Sin{7 t] 11346613 Sin(9 t] 


+ Ce EEE EEE Een a 
206438400 589824000 371589120 


883481 Sm[l1t] 89371 ee 
103219200 137625600 
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Upper = Cos[3t] + Sm[3 t] + 











3Cos[t] 3 3 3Smf[t} 359 | eS o{ 329Cosi|t] 
({- + — Cos[3t] + —— Cos[St] + + —— Sm[3 t] - — Sint 5t]} + [--——» 
16 16 240 16 640 

233 Cos[3t] 479Cos[5t] 27 51Sm[t}] 73913 Sin[3 t] I~ 2a 

+ ——___——— —- — Cos[7t] + + + = Sin'5t} - — Sint 7} + 
1280 1280 640 128 57600 256 640 

3( 89813 Cos[t] 76151 Cos[3t] 139 67 Cos|7 t] 
| s[5t] + —————. - 

153600 230400 600 2048 

469Cos[9t] 831 Smt]  82393Sm[3t] 679Sm[St]  5311Sm[7t]  469Sm[9t] 

ee ee -O 

46080 2048 57600 5120 51200 46080 


‘{ 385217 Cos[t] 105987247 Cos[3 t] 


+ 
614400 309657600 
389277 Cos[5t] 30401 Cos[7t] 327931 Cos[9t] 547Cos[{1]t] 273059Sm{t] 


oe 
1638400 409600 11059200 . 215040 614400 

37616912581 Sm[3t}]  18809Sm[5t] 249571 Sm[7t]  27691Sm[9t] 547 Snf11t] 
23224320000 327680 6144000 2211840 215040 


5 ( 42447078481 Cos[t] 33657149027 Cos[3t] 9222921259 Cos[5t] 1330129 Cos[7 t] 


foe + 
61931520000 92897280000 7 30965760000 39321600 
46069577 Cos[9t] 88301 Cos[]1t] 89371 Cos[13t] 400303339 Sm[t] 82024682603 Sm[3 t] 


ht tt Ht 
9289728000 20643840 137625600 825753600 46448640000 
13130773 Sm{5 t] 16057111 Sm{7 t] 11346613 Sm[9t] 883481 Smf11 t] 89371 Sm 13 t] 


206438400 589824000 371589120 7 103219200 137625600 





APPENDIX E: STABILITY TRANSITION CURVES 


This appendix includes the stability transition curves found in Chapter 5. Curves are given 
analytically through order 10 and numerically through order 30. 


Curves out of A,=0: 


e @€& 8 Te 365 247 <6 11657 <7 1224811 <8 
LowercA = -— +—4+---—;+-—_—_—_—- . ron Oe ee 
2 8&8 32 384 #11024 442368 21233664 2548039680 
65987129 <9 248740367 «10 


611529523200  36691771392000 


Lower sA = 
~0.5¢€ + 0.125¢7 + 0.03125? — 0.018229 «4 + 0.0029297 €° — 
0.00055836 c° — 0.00054899 «7 + 0.00048069 <8 — 0.00010791 €? + 6.7792 x 107° cl? + 
0,000023325 €!! — 0,000020161 «!? + 3.8825 x 107° c!3 + 7.1301 x 1077 «!4 — 1.1828 x 107% el + 
8.9935x 1077 «!®— 1.1027 1077«!” — 9.9766 x 1078 €!8 + 6.0926 x 1078 €!? — 3.6467 x 1078 <9 + 
1.5256x 107! <4! 4 8.7439 1077 7? — 2.8529x 107? € + 9.8924 x 107194 + 3.551 x 107! <# - 
6.2361 x 107! <% + 9.9193 x 107! Fe?” + 2.5033 x 107! 28 — 3.8338 x 1071! <9 + 3.6974 1071! 30 


«e @ 28 Tet 365 247 6 11657¢7 1224811 <8 
UppereA = 4+ > - = - = = Saat Sa, + See 

2 8 32 384 (1024 442368 21233664 2548039680 
65987129 <9 248740367 10 


611529523200 * 3669177139200 








Upper cA = 
0.5€ + 0.125 €2 — 0.03125 & — 0.018229 «4 — 0.0029297 2 - 
0.00055836 «© + 0.00054899 «7 + 0.00048069 <8 + 0.00010791 €? + 6.7792 x 1076 <1 — 
0.000023325 «!! — 0.000020161 «!2 — 3.8825x 10-6 €!3 + 7.1301 x 1077 €!4 + 1.1828 10-6 €l5 + 
8.9935 x 1077 €!© + 1.1027 x 10-7e!7 — 9.9766 x 1078 €!8 — 6.0926 x 10-8 €!9 — 3.6467 x 10-8 «29 — 
1.5256 x 1071! <2! + 8.7439x 10-9 «7 + 2.8529x 107? 3 + 9.8924x 10719 <4 - 3.551 x 10719 <5 - 
6.2361 x 107 !% 2 — 9.9193 x 107! 27 + 2.5033 x 1071! e278 + 3.8338 x 107! 2? + 3.6974 x 1071! 30 


Curves out of Ag=3/2: 


2 1094 34196 297305 8 7797193367 € 10 
LowereA = = — te > ooeoeooooaiel oO 
To 3456 1990656 = 573308928 = 82556485632000 


Lower eA = 





66 
1.5— 0.16667 «2 + 0.03 15394 — 0.0017175¢° — 0.00051858 «® + 0.000094447 «10 + 


0.000013 172 «!2 — 5.7415 x 10-© e!4 — 1.1848 x 1077 €!© + 3.354 1077 €!8 — 2.8733 x 1078 29 
1.7574x 107822 + 3.7778 x 107% «4 + 7.2633 x 107 19 6% — 3.3523 x 107!% 28 — 1.1097 x 107!! ©30 


2 53«4 = 6 983. ® 740213 <8 6745666577 «!0 


— —_— 


e 
=a ee eo 
3 3456 1990656 i 2866544640 = 82556485632000 








3 
UppereA = 5 + 


Upper eA = 
1.54 0.33333 €2 — 0.015336€4 — 0.0035079 © + 0.00025822 €8 + 0.00008171 «1° — 
0.000015556 ¢!2 — 1.5903 x 107° €!4 + 1.0185 10-%e!® — 3.4918 x 1078 €!8 — 5.8479 x 107-8 <2? + 
8.9289 x 1079 22 + 2.7575x 1079 4 ~ 9.066 x 10719 «26 — 8.0056x 10711 <28 + 7.1452 x 107}! _30 


Curves out of Ag=4: 


9 52 WWiTe* 2537 54707€° 11663959 «7 
LowereA = 4- ——--—— + $e > 
a2 64 40960 = 737280) =: 15728640 ~—-:14155776000 
17260922723 «8 43486217423 €? 2804 162996941 «10 


12683575296000  1223059046400000 . 14611478740992000 





Lower sA = 
4. — 0.28125¢* — 0.078125e3 + 0.070239 «4 + 0.003441? — 0.0034782«° + 
0,00082397 «” — 0.0013609«® — 0,.000035555€” + 0,00019192«!° — 0,000016118¢!! + 
0.000037433 «!? + 1.4085 x 107° «3 — 0,000011912 «!* + 2.053 x 1077! — 6.2434 1077 €!6 — 
6.8113 x 10-8 <!” + 7.0367 107718 + 1.5598 x 1078 !? — 4.6186 x 1078” + 2.0464 107? «2! — 
3.6947 x 107-8 <2 — 2.1736 1077 e+ 7.4707 x 107? «4 + 9.323 x 1071 + 1.5107 107? e* + 
1.8088 x 107! <2? ~ 6.9398 x 107! <?8 - 2.4687 x 1071! <? — 2.0128 x 107!! <0 


9e2 553) 2877e4 = 253765 =: 54707 €® 11663959 «7 
UppereA = 4- —— + — + ——— 
32 64 40960 737280 15728640 14155776000 

17260922723 «8 43486217423 &? 2804 16299694 1 «19 


12683575296000  1223059046400000  14611478740992000 


RR em eee 


Upper eA = 

4,— 0.28125? + 0.078125 3 + 0.070239 €4 — 0.003441 €> — 0.0034782«6 — 
0,00082397 «7 — 0.0013609 «8 + 0.000035555 €? + 0,00019192 «!9 + 0.000016118«!! + 
0.000037433 «!2 — 1.4085 x 107° €}3 — 0,000011912 €!4 — 2.053 x 107-715 - 6.2434 1077 «16 + 
6.8113 x 1078!” + 7.0367 1077 €!8 — 1.5598x 1078 €!9 — 4.6186 x 10-8 <29 — 2.0464 1079 «2! - 
3.6947 x 1078 «22 + 2.1736 1077 23 + 7.4707 x 1079 e*4 — 9.323 x 107! 25 + 1.5107 1079 ® - 
1.8088 x 1071 €27 — 6.9398 x 1071928 + 2.4687 x 107}! 29 — 2.0128 x 1071! <30 


on SineS 4 qe. 
=a oust. 
1 199225277187 3 
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